!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! DTM FOR 3D PROBLEM
!!-------------------
!! written by Vivien R. Lecoustre (vlecous1@umd.edu)
!! This file contains three modules:
!!    radiation_m 
!!    dtm_m
!!    topology_m
!! To use, include this file in your makefile.
!! Make sure your compiler is compatible Fortran 2003
!----------------------------------------------------------------------------------------
!
! How to use:
!------------
! !!This package requires the user to provide information from an external module named:
!       param_m.
! The information needed from param_m are:
! nx,ny,nz :: Type Integer. 
!             Corresponds to the number of mesh points in x, y, and z direction
!  x, y, z :: type Real. Dimension: x(nx), y(ny), z(nz). These are the mesh point location
!             in dimensionless unit
!  l_ref   :: Real. Coefficient to used to make x, y, z. non dimensional
! periodic_x, periodic_y, periodic_z :: Type Integer. 
!             if 0, non-periodic BC. If 1, periodic BC for radiation
!
! User input needed.
! The user has to provide input into a file named radiation.in
! Structure of radiation.in:
! i_radiation
! i_sad     method discretization solid angle. 0 = constant discrete solid angle. if 0 provides n_rays
! n_rays        
! n_phi
! n_theta
! emiss_t
! emiss_b
! emiss_l 
! emiss_r
! emiss_s
! emiss_n
!
! II Sequence for use:
!---------------------
! 1) call initialize_radiation
! 2) call DTM
!   After 2 the value of the divergence for radiative source term is calculated.
!   radiation_m provides a list of post-processing subroutines that will write tecplot files
!   |Print_irrad
!   |Print_Fluxes
!   |Print_Cell_Source
!   |Print_Diag_North
!   |Print_Div_qr(div_qr)
! 3) when done, call allocate_radiation_arrays(-1) the flag -1 is very important.
!    This will free some memory.
!----------------------------------------------------------------------------------------

!=========================================================================================
  module radiation_m
!=========================================================================================
! module for Discrete Transfer Method. This version is for validation test only.
! Do not use for Parallel work.
! Works for non-uniform grid 
!----------------------------------------------------------------------------------------
! Based on the Yi Wang version (04/20/04)
! Yi Wang (yiwang@wam.umd.edu)
!----------------------------------------------------------------------------------------
! Adapted by Vivien Lecoustre (vlecous1@umd.edu)
! 03/12/12 
!=========================================================================================
! Transformation Matrices P_ij - USED IN DTM SUBROUTINE
! (south) -> (north)  (x,y,z) -> (x,-y,-z)
!       |1  0  0|   
! P_SN =|0 -1  0|  
!       |0  0 -1| 

! (south) -> (right)  (x,y,z) -> (y,z,x)
!       |0  1  0|   
! P_SR =|0  0  1|  
!       |1  0  0| 

! (south) -> (left) (x,y,z) -> (y,-z,-x)
!       |0 -1  0|   
! P_SL =|0  0  1|  
!       |-1 0  0| 

! (south) -> (bottom) (x,y,z) -> (z,x,y)
!       |0  0  1|   
! P_SB =|1  0  0|  
!       |0  1  0| 

! (south) -> (top)  (x,y,z) -> (z,-x,-y)
!       |0  0  1|   
! P_ST =|-1 0  0|  
!       |0 -1  0| 
!-----------------------------------From North Plane-------------------------------------
! (north) -> (south)  (x,-y,-z) -> (x,y,z)
!       |1  0  0|   
! P_NS =|0 -1  0| (=P_SN) 
!       |0  0 -1| 

! (north) -> (right)  (x,-y,-z) -> (y,z,x)
!       |0 -1  0|   
! P_NR =|0  0 -1|  
!       |1  0  0| 

! (north) -> (left) (x,-y,-z) -> (-y,z,-x)
!       |0  1  0|   
! P_NL =|0  0 -1|   
!       |-1 0  0| 

! (north) -> (bottom) (x,-y,-z) -> (z,x,y)
!       |0  0 -1|   
! P_NB =|1  0  0|  
!       |0 -1  0| 

! (north) -> (top)  (x,-y,-z) -> (z,-x,-y)
!       |0  0 -1|   
! P_NT =|-1 0  0|  
!       |0  1  0| 
!--------------------------From Right Plane----------------------------------------------
! (right) -> (south)  (y,z,x) -> (x,y,z)
!       |0  0  1|   
! P_RS =|1  0  0| (=P_SB)  
!       |0  1  0| 

! (right) -> (north)  (y,z,x) -> (x,-y,-z)
!       |0  0  1|   
! P_RN =|-1 0  0| (=P_ST) 
!       |0 -1  0| 

! (right) -> (left) (y,z,x) -> (-y,z,-x)
!       |-1 0  0|   
! P_RL =|0  1  0|   
!       |0  0 -1| 

! (right) -> (bottom) (y,z,x) -> (z,x,y)
!       |0  1  0|   
! P_RB =|0  0  1| (=P_SR)
!       |1  0  0| 

! (right) -> (top)  (y,z,x) -> (z,-x,-y)
!       |0  1  0|   
! P_RT =|0  0 -1| (=P_NL)  
!       |-1 0  0| 

!--------------------------From Left Plane----------------------------------------------
! (left) -> (south)  (-y,z,-x) -> (x,y,z)
!       |0  0  1|   
! P_LS =|-1 0  0|  (=P_ST)
!       |0  -1 0| 

! (left) -> (north)  (-y,z,-x) -> (x,-y,-z)
!       |0  0 -1|   
! P_LN =|1  0  0| (=P_NB)
!       |0 -1  0| 

! (left) -> (right) (-y,z,-x) -> (y,z,x)
!       |-1 0  0|   
! P_LR =|0  1  0|  (=P_RL)
!       |0  0 -1| 

! (left) -> (bottom) (-y,z,-x) -> (z,x,y)
!       |0  1  0|   
! P_LB =|0  0 -1| (=P_NL)
!       |-1 0  0| 

! (left) -> (top)  (-y,z,-x) -> (z,-x,-y)
!       |0  1  0|   
! P_LT =|0  0  1| (=P_SR)  
!       |1  0  0| 

!--------------------------From Bottom Plane----------------------------------------------
! (bottom) -> (south)  (z,x,y) -> (x,y,z)
!       |0  1  0|   
! P_BS =|0  0  1|  (=P_SR)   
!       |1  0  0| 

! (bottom) -> (north)  (z,x,y) -> (x,-y,-z)
!       |0  1  0|   
! P_BN =|0  0 -1|  (=P_NL)
!       |-1 0  0| 

! (bottom) -> (right) (z,x,y) -> (y,z,x)
!       |0  0  1|   
! P_BR =|1  0  0|  (=P_SB)
!       |0  1  0| 

! (bottom) -> (left) (z,x,y) -> (-y,z,-x)
!       |0  0 -1|   
! P_BL =|1  0  0| (=P_NB)
!       |0 -1  0| 

! (bottom) -> (top)  (z,x,y) -> (z,-x,-y)
!       |1  0  0|   
! P_LT =|0 -1  0| (=P_SN)  
!       |0  0 -1| 

!--------------------------From Top Plane----------------------------------------------
! (top) -> (south)  (z,-x,-y) -> (x,y,z)
!       |0 -1  0|   
! P_TS =|0  0 -1|  (=P_NR)   
!       |1  0  0| 

! (top) -> (north)  (z,-x,-y) -> (x,-y,-z)
!       |0 -1  0|   
! P_TN =|0  0  1| (=P_SL)
!       |-1 0  0| 

! (top) -> (right) (z,-x,-y) -> (y,z,x)
!       |0  0 -1|   
! P_TR =|1  0  0|  (=P_NB)
!       |0 -1  0| 

! (top) -> (left) (z,-x,-y) -> (-y,z,-x)
!       |0  0  1|   
! P_TL =|1  0  0|  (=P_SB) 
!       |0  1  0| 

! (top) -> (bottom)  (z,-x,-y) -> (z,x,y)
!       |1  0  0|   
! P_TB =|0 -1  0| (=P_SN)  
!       |0  0 -1| 
!---------------------------------------------------------------------------------------- 
!
!  Recall: 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top
!
  implicit none

!-----------------------------------------------------------------------------------------
! type definition
  
!  private

!   type path
! !    integer,dimension(:),pointer :: ix         !x index of control volume that a ray passed
! !    integer,dimension(:),pointer :: iy         !y index of control volume that a ray passed
! !    integer,dimension(:),pointer :: iz         !z index of control volume that a ray passed
! 
!     type(cuboid), dimension(:), pointer :: Passed_CV  ! New. Will replace ix, iy, and iz to get adresses of controlled
!                                         ! volume crossed
! 
!     real,  dimension(:),pointer :: length     !beam length of ray in the control volume
! !    real   ,dimension(:),pointer :: Surface   !Area of the control volume crossed by the ray
! !   Obsolete, use cell.Surface instead. Saves some memory
!   end type path

!----------------------------------------------------------------------------------------
! Version control variables
  character(255), parameter, private :: Radid   ='$Id: radiation_m.f90 34 2013-07-23 21:43:05Z vlecous1 $'
  character(255), parameter, private :: Radrev  ='$Revision: 34 $'
  character(255), parameter, private :: Raddate ='$Date: 2013-07-23 21:43:05 +0000 (Tue, 23 Jul 2013) $'

!----------------------------------------------------------------------------------------
! Define Fortran derived types

  type cuboid_pntr
    type(cuboid), pointer :: cub
  end type cuboid_pntr

  type coordinates
    real :: x
    real :: y
    real :: z
  end type coordinates

  type facet
    type(coordinates), dimension(4) :: vertex
    type(coordinates)               :: center
    type(cuboid), pointer           :: cell
    real                            :: Surface
    real                            :: Temp      ! Temperature of the facet
    integer                         :: i_dim1    ! Indices of the first dimension
    integer                         :: i_dim2    ! Indices of the second dimension
    type(face), pointer             :: boundary  ! pointer to the boundary face if facet belongs to a boundary
  end type facet 
    
  type face
    type(facet), allocatable, dimension(:,:) :: element
    integer                                  :: fixed_dim
    real                                     :: Fixed_coordinate
    real                                     :: emissivity
    integer                                  :: plane
    integer                                  :: coordinate
    type(coordinates), dimension(4)          :: vertex ! location of the verteces of the face
  end type face

  type cuboid
    real, dimension(2) :: x                    ! x coordinates of the cuboid
    real, dimension(2) :: y                    ! y coordinates of the cuboid
    real, dimension(2) :: z                    ! z coordinates of the cuboid
    real               :: Surface              ! Surface of the control volume
    real               :: Volume
    real               :: Temp                 ! Temperature of the control volume (in K)
    real               :: c_abs                ! Absorption coefficient for the control volume
    real               :: steradiancy          ! BlackBody radiative intensity per unit of steradiant
    real               :: Source_cv            ! Radiative source for the control volume
    integer            :: i_x                  ! Indice of the cell in x direction
    integer            :: i_y                  ! Indice of the cell in y direction
    integer            :: i_z                  ! Indice of the cell in z direction
  end type cuboid

  type beam
    type(facet), pointer                         :: Facet_origin        
    type(cuboid_pntr), dimension(:), allocatable :: Passed_CV  ! get adresses of control volume crossed
! For debugging only. REMOVE AFTER
    type(coordinates), dimension(:), allocatable :: Point      ! Point crossed by ray_buffer ! TO REMOVE IN FINAL VERSION
! END
    real, allocatable, dimension(:)              :: Length     ! beam length of ray in the control volume
! For debugging only. REMOVE AFTER
    real    :: Temp_origin
! END 
    integer :: n_path    
    integer :: iray2
    integer :: iphi
    integer :: itheta
    integer :: edge_final
  end type beam

  type Radiance_Boundary
    real, allocatable :: South(:,:,:)
    real, allocatable :: North(:,:,:)
    real, allocatable :: Right(:,:,:)
    real, allocatable :: Left(:,:,:)
    real, allocatable :: Bottom(:,:,:)
    real, allocatable :: Top(:,:,:)
  end type Radiance_Boundary

  type Angle_Azimuthal
    real, allocatable :: azimuth(:)
  end type Angle_Azimuthal

  type(Radiance_Boundary),target :: rib             ! Variable containing the incident intensities
                                                    ! on each boundary points  

  type(face), target, dimension(6) :: Boundaries        ! Array of object faces. Contains properties of the 6 domain faces

  type(beam),   allocatable, target  :: ray(:,:)        ! order of beams passing CVs along each ray
  type(cuboid), allocatable, target  :: cell_tot(:,:,:) ! Array containig all the control volumes
                                                          ! dimensions should be mx-1, my-1, mz-1

  type(cuboid) :: Global_Cell ! The computational domain is defined as a cuboid
                              ! Variable needed for functions: findpath3d_new, Get_Facet_Index

!----------------------------------------------------------------------------------------
! integers

  integer :: i_radiation ! 0-no radiation, 1-DTM, 2-DOM, 3-optically thin 
  integer :: mx, my, mz  ! number of grid points including ghost cells (nx+1 or nx)
  integer :: n_points    ! number of boundary cells having incident rays 
                         ! (2*(mx-1)*(my-1)+2*(my-1)*(mz-1)+2*(mx-1)*(mz-1))
! Solid angle discretization
  integer :: i_sad       !Method discretization of solid angles. 
                         ! if 0 = constant discrete solid angles. User provides n_rays
                         ! if 1 = constant delta_theta and delta_phi

  integer :: n_theta     !Decomposition of theta, polar angle. Number of rays fired by
                         !each boundary point at a given angle phi.     

  integer, allocatable :: n_phi(:) ! Decomposition of phi, azimuthal angle. Number of rays fired by
                                   ! each boundary point for a given angle tetha. Allow non uniform
                                   ! decomposition of azimuthal angle
  
  type(Angle_Azimuthal), allocatable, target :: phi(:)     ! phi: Azimutal angle, spans from 0 to 2 pi (excluded)
  type(Angle_Azimuthal), allocatable, target :: cos_int(:) ! integral of the cosine of angle between ray
                                                   ! and surface normal over element solid angle

  real, allocatable :: theta(:)                  ! theta: Polar angle, spans from 0 to pi/2 (included)
                         
  integer :: n_rays 	 ! Number of ray fired by each boundary points,
                         ! if i_sad = 1 => n_rays = n_phi*n_theta
                         ! if i_sad = 0 => n_rays = sum(n_phi)

  integer :: output_unit ! fortran file unit for output message. Needed for error message

!----------------------------------------------------------------------------------------
! reals  

  real,allocatable :: rad_source(:,:,:)  !radiative source at grid points
  real,allocatable,target :: xx(:)       !coordinate including ghost cell (in SI unit)
  real,allocatable,target :: yy(:)       !coordinate including ghost cell (in SI unit)
  real,allocatable,target :: zz(:)       !coordinate including ghost cell (in Si unit)

  real,allocatable :: irrad0(:)          !irradiation on the boundary points (old value)
  real,allocatable :: irrad(:)           !irradiation on the boundary points (new value)
  real,allocatable :: emission(:)        !emissive power on the boundary points

  real,allocatable :: kappa(:,:,:)       !absorption coefficient at the grid point

  real,allocatable,target :: rivb_t(:,:,:) !Radiative Intensity at the Virtual Boundary, y=ymax
  real,allocatable,target :: rivb_b(:,:,:) !Radiative Intensity at the Virtual Boundary, y=ymin
  real,allocatable,target :: rivb_l(:,:,:) !Radiative Intensity at the Virtual Boundary, x=xmax
  real,allocatable,target :: rivb_r(:,:,:) !Radiative Intensity at the Virtual Boundary, x=xmin
  real,allocatable,target :: rivb_s(:,:,:) !Radiative Intensity at the Virtual Boundary, z=zmin
  real,allocatable,target :: rivb_n(:,:,:) !Radiative Intensity at the Virtual Boundary, z=zmax

  real,allocatable,target :: rivb0_t(:,:,:) !Radiative Intensity at the Virtual Boundary, y=ymax
  real,allocatable,target :: rivb0_b(:,:,:) !Radiative Intensity at the Virtual Boundary, y=ymin
  real,allocatable,target :: rivb0_l(:,:,:) !Radiative Intensity at the Virtual Boundary, x=xmax
  real,allocatable,target :: rivb0_r(:,:,:) !Radiative Intensity at the Virtual Boundary, x=xmin
  real,allocatable,target :: rivb0_s(:,:,:) !Radiative Intensity at the Virtual Boundary, z=zmin
  real,allocatable,target :: rivb0_n(:,:,:) !Radiative Intensity at the Virtual Boundary, z=zmax

  real,parameter :: sigma = 5.670373e-8  ![W/(m^2 * K^4)] = [kg/(s^3 * K^4)] (NIST Value, 2011)

  real,parameter :: criterion = epsilon(sigma) ! Criterion used in geometric functions to assess
                                               ! condition of belonging in cell
  real,parameter :: criterion_face = criterion*1.0d+4 ! Criterion used in geometric functions to assess
                                                      ! condition of belonging in facet

  real :: emiss_t, emiss_b                ! emissivity of y=ymax plane, emissivity of y=ymin plane
  real :: emiss_l, emiss_r                ! emissivity of x=xmax plane, emissivity of x=xmin plane
  real :: emiss_s, emiss_n                ! emissivity of z=zmin plane, emissivity of z=zmax plane

  real :: Temp_South   ! Fixed temperature for South boundary, in Kelvin
  real :: Temp_North   ! Fixed temperature for North boundary, in Kelvin
  real :: Temp_Right   ! Fixed temperature for South boundary, in Kelvin
  real :: Temp_Left    ! Fixed temperature for South boundary, in Kelvin
  real :: Temp_Bottom  ! Fixed temperature for South boundary, in Kelvin
  real :: Temp_Top     ! Fixed temperature for South boundary, in Kelvin

! Variable tensor Change_plane(6,6,3,3). Used in DTM subroutine to compute change of basis needed
! to fetch the radiation intensities at the virtual bondary
  
  real, dimension(6,6,3,3) :: Change_plane

  real, parameter :: pi = acos(-1.0)

!----------------------------------------------------------------------------------------
! Variables for MPI Cartesian coordinates. Assume 3D decomposition
  logical, dimension(3) :: periodic_bc ! For Cartesian domain decomposition
  integer, dimension(3) :: cart_dim    ! # of processes in each dimension
  integer, dimension(3) :: myid_cart   ! Coordinates of the current process in cart_comm

  integer  :: cart_comm        ! MPI Cartesian topology communicator
  integer  :: size_type_real   ! Size of real number (use for MPI communication)
  integer  :: TypeReal         ! MPI DATATYPE representing a real number
  integer  :: ierr_MPI_rad

!----------------------------------------------------------------------------------------
! Logical for need of ghost plan
  logical :: ghost_x, ghost_y, ghost_z

! Logical array for Virtual boundaries
  logical, dimension(6) :: is_vb   ! TRUE = virtual boundary, FALSE = real boundary

!----------------------------------------------------------------------------------------
! Overloading of some procedures

  interface metric
    module procedure metric_1d
    module procedure metric_3d
  end interface

  interface surface
    module procedure surface_2D, surface_3D
  end interface

!----------------------------------------------------------------------------------------
! Subroutines and variables attributes:
! Public (no restriction) or Private (restricted to within this module)

  public  :: initialize_radiation, allocate_radiation_arrays, Update_Cell,              &
             Get_Points_Index, Transf_Ptwise_CV, Radiative_emission, exchange_RIVB,     &
             Get_Div_qr, tecplot_write_xy_line, print_ray, Print_Diag_North,            &
             Print_Cell_Source, Print_Div_qr, Print_Fluxes, Print_irrad, Init_MPI_Rad,  &
             metric, Print_all_rays

  private :: Error_MSG, findpath3d_new, init_rad_geometry, initialize_ChangePlane,      &
             create_CV, get_absorption, Find_Associated_Ray, metric_1d, metric_3d

  public  :: pi, rad_source, Change_plane, sigma

  private :: i_sad

!-----------------------------------------------------------------------------------------
  contains
!=========================================================================================
  subroutine initialize_radiation(io)
!=========================================================================================
  use param_m,    only : numdim,nx,ny,nz, x, y, z, l_ref
  use topology_m
  
  implicit none

  integer, intent(in) :: io

  integer :: i, j
  integer :: n_phi1
  real :: del_theta, del_phi

! MPI Local Variables

  integer, dimension(3) :: coordinates
  integer  :: target_id_p, target_id_m ! Rank of remote processed. _m (minus) _p (plus)

  integer  :: req          ! Needed by MPI
  integer  :: my_px_id, my_py_id, my_pz_id

  character(LEN=100):: filename

! Assign value to output_unit
  output_unit = io
!----------------------------------------------------------------------------------------
! write header

  if (myid==0) then
    filename='radiation.in'
    open(unit=1,file=trim(filename),status='old')
    read(1,*) i_radiation
    close (1)
  endif
 
  call MPI_Bcast(i_radiation, 1, MPI_INTEGER, 0, gcomm, ierr_MPI_rad)

  select case(i_radiation)
    case(0)
      if (myid.eq.0) then
        write(io,*) 'Radiation module is turned off'
        write(io,*)
      endif
      return ! Exit radiation module
    case(1)
      if (myid.eq.0) then
        write(io,*) 'Initializing radiation module: 3D Discrete Transfer Method'
        write(io,*)
      endif
    case default
      if (myid.eq.0) then
        write(io,*) 'Wrong flag for radiation module'
        write(io,*)
      endif
      return ! Exit radiation module
  end select

! consistancy check

  if (numdim.ne.3) then
    if(myid.eq.0) then
      write(io,*) 'current DTM radiation module only works for 3D problems'
      write(io,*)
    endif    
  endif

!----------------------------------------------------------------------------------------   
! Read the fitting coefficients for Planck mean absorption coefficients 
! borrow the routine in opthinrad_m

  call initialize_ChangePlane()

!----------------------------------------------------------------------------------------
! Initialize parameters mx, my
! note: mx, my, mz are local point numbers for radiation (including ghost cell) 

  if (ghost_x) then 
    mx = nx+1
  else
    mx = nx
  endif

  if (ghost_y) then 
    my = ny+1
  else
    my = ny
  endif

  if (ghost_z) then 
    mz = nz+1
  else
    mz = nz
  endif
  
!----------------------------------------------------------------------------------------
! Read input parameters for DTM

  if (myid==0) then
    open(unit=1,file=trim(filename),status='old')
    read(1,*) 
    read(1,*) i_sad
    read(1,*) n_rays
    read(1,*) n_phi1
    read(1,*) n_theta
    read(1,*) emiss_t
    read(1,*) emiss_b
    read(1,*) emiss_l 
    read(1,*) emiss_r
    read(1,*) emiss_s
    read(1,*) emiss_n
    close (1)
  endif

  call MPI_Bcast(i_sad   , 1, MPI_INTEGER, 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(n_rays  , 1, MPI_INTEGER, 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(n_phi1  , 1, MPI_INTEGER, 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(n_theta , 1, MPI_INTEGER, 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(emiss_t , 1, MPI_REAL8  , 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(emiss_b , 1, MPI_REAL8  , 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(emiss_l , 1, MPI_REAL8  , 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(emiss_r , 1, MPI_REAL8  , 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(emiss_s , 1, MPI_REAL8  , 0, gcomm, ierr_MPI_rad)
  call MPI_Bcast(emiss_n , 1, MPI_REAL8  , 0, gcomm, ierr_MPI_rad)

!----------------------------------------------------------------------------------------

  if (i_sad == 0) then 
    n_theta = ceiling(1.17*float(n_rays)**(1.0/2.26))
! n_theta should be an even number
    n_theta = n_theta - mod(n_theta,2)
  endif

  del_theta= 0.5*pi/real(n_theta) ! theta: Polar angle, spans from 0 to pi/2

  allocate(theta(n_theta))
  allocate(n_phi(n_theta))
  allocate(phi(n_theta))
  allocate(cos_int(n_theta))

  do concurrent (j = 1:n_theta)

    theta(j) = (dfloat(j)-0.5)*del_theta

    if (i_sad == 0) then   
      n_phi(j) = ceiling(max(4.0,float(n_rays)*(sin(theta(j))*sin(del_theta))))
      n_phi(j) = n_phi(j) - mod(n_phi(j),4)
    else
      n_phi(j) = n_phi1
    endif

    del_phi  = 2.0*pi/real(n_phi(j))

    allocate(phi(j)%azimuth(n_phi(j)))
    allocate(cos_int(j)%azimuth(n_phi(j)))

    do i = 1, n_phi(j)
       phi(j)%azimuth(i)     = (dfloat(i)-0.5)*del_phi
       cos_int(j)%azimuth(i) = del_phi*sin(theta(j))*cos(theta(j))*sin(del_theta) 
    enddo
  enddo

  n_rays = sum(n_phi)

  if (myid==0) then
    write(io,*) 'Number of rays: ', n_rays
    write(io,*) 'Number of theta discretization: ', n_theta
    write(io,*) 'Number of phi discretization: ', n_phi
  endif

! Define number of points for intensity vector
  n_points=2*(mx-1)*(my-1)+2*(my-1)*(mz-1)+2*(mx-1)*(mz-1)

  call allocate_radiation_arrays(1)

!----------------------------------------------------------------------------------------
! grid locally used by radiation model

  xx(1:nx)=x(1:nx)*l_ref
  yy(1:ny)=y(1:ny)*l_ref
  zz(1:nz)=z(1:nz)*l_ref

!----------------------------------------------------------------------------------------
! Get ghost point
    
  my_px_id = myid_cart(1)
  my_py_id = myid_cart(2)
  my_pz_id = myid_cart(3)

! Get ghost point in x
!---------------------

  if( my_px_id.gt.0 ) then
! Get address (in target_disp_m) of processor in location (my_px_id-1,my_py_id,my_pz_id)
! within communicator cart_comm

    coordinates(1) = my_px_id-1
    coordinates(2) = my_py_id
    coordinates(3) = my_pz_id
    target_id_m    = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id_m,ierr_MPI_rad)  

! Send the first element of x to target_id_m
    call MPI_Send(x(1)*l_ref,1,MPI_REAL8,target_id_m,1,cart_comm,ierr_MPI_rad)
  endif
 
  if( my_px_id.lt.cart_dim(1)-1 ) then
! Get address (in target_disp_m) of processor in location (my_px_id+1,my_py_id,my_pz_id)
! within communicator cart_comm
 
    coordinates(1) = my_px_id+1
    coordinates(2) = my_py_id
    coordinates(3) = my_pz_id
    target_id_p    = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id_p,ierr_MPI_rad)  

! Send the last element of xx from target_id_p
    call MPI_Recv(xx(mx),1,MPI_REAL8,target_id_p,1,cart_comm,status,ierr_MPI_rad)
  endif

! Get ghost point in y
!---------------------

  if( my_py_id.gt.0 ) then
! Get address (in target_disp_m) of processor in location (my_px_id,my_py_id-1,my_pz_id)
! within communicator cart_comm

    coordinates(1) = my_px_id
    coordinates(2) = my_py_id-1
    coordinates(3) = my_pz_id
    target_id_m    = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id_m,ierr_MPI_rad)  

! Send the first element of y to target_id_m
    call MPI_Send(y(1)*l_ref,1,MPI_REAL8,target_id_m,2,cart_comm,ierr_MPI_rad)
  endif
 
  if( my_py_id.lt.cart_dim(2)-1 ) then
! Get address (in target_disp_m) of processor in location (my_px_id,my_py_id+1,my_pz_id)
! within communicator cart_comm
 
    coordinates(1) = my_px_id
    coordinates(2) = my_py_id+1
    coordinates(3) = my_pz_id
    target_id_p    = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id_p,ierr_MPI_rad)  

! Send the last element of yy from target_id_p
    call MPI_Recv(yy(my),1,MPI_REAL8,target_id_p,2,cart_comm,status,ierr_MPI_rad)
  endif


! Get ghost point in z

!---------------------
  if( my_pz_id.gt.0 ) then
! Get address (in target_disp_m) of processor in location (my_px_id,my_py_id,my_pz_id-1)
! within communicator cart_comm

    coordinates(1) = my_px_id
    coordinates(2) = my_py_id
    coordinates(3) = my_pz_id-1
    target_id_m    = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id_m,ierr_MPI_rad)  

! Send the first element of y to target_id_m
    call MPI_Send(z(1)*l_ref,1,MPI_REAL8,target_id_m,3,cart_comm,ierr_MPI_rad)
  endif
 
  if(my_pz_id.lt.cart_dim(3)-1) then
! Get address (in target_disp_m) of processor in location (my_px_id,my_py_id,my_pz_id+1)
! within communicator cart_comm
 
    coordinates(1) = my_px_id
    coordinates(2) = my_py_id
    coordinates(3) = my_pz_id+1
    target_id_p    = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id_p,ierr_MPI_rad)  

! Send the last element of yy from target_id_p
    call MPI_Recv(zz(mz),1,MPI_REAL8,target_id_p,3,cart_comm,status,ierr_MPI_rad)
  endif

!----------------------------------------------------------------------------------------
!                           Create control volume arrays
  call create_CV()

!----------------------------------------------------------------------------------------
!                           Create Boundaries Surfaces
  call Create_facets()

!----------------------------------------------------------------------------------------
! get geometry information for rays

  call init_rad_geometry(mx,my,mz,xx,yy,zz,io)

!----------------------------------------------------------------------------------------
! write header

  if(myid.eq.0) then
    write(io,*) 'Radiation module for DTM initialization finished.'    
  endif
!----------------------------------------------------------------------------------------
  return
  end subroutine initialize_radiation

!=========================================================================================
  subroutine allocate_radiation_arrays(flag)
!=========================================================================================
! allocate radiation arrays
! Initialize HDF5 for flag=1
! Finalize HDF5 and MPI when flag =-1
!-----------------------------------------------------------------------------------------
  use param_m,    only : nx, ny, nz
  use topology_m, only : ierr, hdf5_error
  use HDF5
  implicit none
!-----------------------------------------------------------------------------------------
! declarations passed

  integer flag
!-----------------------------------------------------------------------------------------
! radiation arrays

  if(flag.eq.1) then

    allocate(rad_source(nx,ny,nz)); rad_source = 0.0
    allocate(ray(n_rays,n_points))         !note: n_points is at the outer loop

    allocate(irrad0(n_points))
    allocate(irrad(n_points))
    allocate(emission(n_points))

!--------Allocate rib variables------------------------------------------------------------
! Recall:
! South   & North are mx-1 by my-1 points
! Left   & Right  are my-1 by mz-1 points
! Bottom & Top    are mx-1 by my-1 points

    allocate(rib%South( n_rays,mx-1,my-1))  
    allocate(rib%North( n_rays,mx-1,my-1))  
    allocate(rib%Left(  n_rays,my-1,mz-1))   
    allocate(rib%Right( n_rays,my-1,mz-1))  
    allocate(rib%Bottom(n_rays,mx-1,mz-1)) 
    allocate(rib%Top(   n_rays,mx-1,mz-1))    
    
    rib%South  = 0.0
    rib%North  = 0.0
    rib%Right  = 0.0
    rib%Left   = 0.0
    rib%Bottom = 0.0
    rib%Top    = 0.0

    allocate(rivb_t(n_rays,mx-1,mz-1)) ; rivb_t=0.0
    allocate(rivb_b(n_rays,mx-1,mz-1)) ; rivb_b=0.0
    allocate(rivb_l(n_rays,my-1,mz-1)) ; rivb_l=0.0
    allocate(rivb_r(n_rays,my-1,mz-1)) ; rivb_r=0.0
    allocate(rivb_n(n_rays,mx-1,my-1)) ; rivb_n=0.0
    allocate(rivb_s(n_rays,mx-1,my-1)) ; rivb_s=0.0

    allocate(rivb0_t(n_rays,mx-1,mz-1)) ; rivb0_t=0.0
    allocate(rivb0_b(n_rays,mx-1,mz-1)) ; rivb0_b=0.0
    allocate(rivb0_l(n_rays,my-1,mz-1)) ; rivb0_l=0.0
    allocate(rivb0_r(n_rays,my-1,mz-1)) ; rivb0_r=0.0
    allocate(rivb0_n(n_rays,mx-1,my-1)) ; rivb0_n=0.0
    allocate(rivb0_s(n_rays,mx-1,my-1)) ; rivb0_s=0.0

    allocate(xx(mx))
    allocate(yy(my))
    allocate(zz(mz))

    allocate(cell_tot(mx-1,my-1,mz-1)) 

    allocate(kappa(nx,ny,nz)); kappa=0.0

! Initialize HDF5
    call h5open_f(hdf5_error)

! Call Init_MPI_Rad
  elseif(flag.eq.-1) then

!----------------------------------------------------------------------------------------
! cleanup MPI processes

    call MPI_Finalize(ierr)

!----------------------------------------------------------------------------------------
! Close HDF5
    call h5close_f(hdf5_error)

    deallocate(rad_source)
    deallocate(ray)
    deallocate(n_phi)
    deallocate(cos_int)
    deallocate(phi)
    deallocate(theta)
    deallocate(irrad0)
    deallocate(irrad)
    deallocate(emission)

    deallocate(rib%South)
    deallocate(rib%North)
    deallocate(rib%Right)
    deallocate(rib%Left)
    deallocate(rib%Bottom)
    deallocate(rib%Top)

    deallocate(rivb_t)
    deallocate(rivb_b)
    deallocate(rivb_l)
    deallocate(rivb_r)
    deallocate(rivb_n)
    deallocate(rivb_s)
    deallocate(xx)
    deallocate(yy)
    deallocate(zz)

    if (allocated(cell_tot))              deallocate(cell_tot)

    if (allocated(Boundaries(1)%element)) deallocate(Boundaries(1)%element)
    if (allocated(Boundaries(2)%element)) deallocate(Boundaries(2)%element)
    if (allocated(Boundaries(3)%element)) deallocate(Boundaries(3)%element)
    if (allocated(Boundaries(4)%element)) deallocate(Boundaries(4)%element)
    if (allocated(Boundaries(5)%element)) deallocate(Boundaries(5)%element)
    if (allocated(Boundaries(6)%element)) deallocate(Boundaries(6)%element)

    if (allocated(kappa))                 deallocate(kappa)
    
  endif
!-----------------------------------------------------------------------------------------
  return
  end subroutine allocate_radiation_arrays
 
!=========================================================================================
  pure function Is_on_cuboid_face(Point,cell)
!=========================================================================================
! This function tests if the point P (type coordinates) belongs to the surface on the cuboid "cell".
! Point belongs to the surface of cell is at least one of the condition are met:
! Test if Point is on one of the face.
! If point belongs to cuboid, than the following statements are true:
! d(cell%x(1),cell%x(2)) = d(x,cell%x(1))+d(x,cell%x(2))
! d(cell%y(1),cell%y(2)) = d(y,cell%y(1))+d(y,cell%y(2))
! d(cell%z(1),cell%z(2)) = d(z,cell%z(1))+d(z,cell%z(2))
! Point belongs to one of the face of the cuboid if at least one
! of its coordinates verifies the hyperplane equation.
! Example, If Point belongs to cuboid cell and has y=y1, then we can affirm that
! point belong to the face of equation y=y1 of the cuboid considered.
! Hence we will test that at least one of the d(cell%i,Point(i)) is less than criterion
! Returns false if not
! The criterion to be on the face is:
! Criterion_face
!----------------------------------------------------------------------------------------

  implicit none

  type(cuboid), intent(in)        :: cell
  type(coordinates), intent(in)   :: Point

  logical :: Is_on_cuboid_face

  Is_on_cuboid_face = .false.
! If P is on the surface of the cuboid, it is also inside the cuboid
  
  if (Is_inside_cuboid(Point,cell)) then 
! First dimension
    if (metric(Point%x,cell%x(1))<= criterion_face) Is_on_cuboid_face = .true.
    if (metric(Point%x,cell%x(2))<= criterion_face) Is_on_cuboid_face = .true.
! Second dimension
    if (metric(Point%y,cell%y(1))<= criterion_face) Is_on_cuboid_face = .true.
    if (metric(Point%y,cell%y(2))<= criterion_face) Is_on_cuboid_face = .true.
! Third dimension
    if (metric(Point%z,cell%z(1))<= criterion_face) Is_on_cuboid_face = .true.
    if (metric(Point%z,cell%z(2))<= criterion_face) Is_on_cuboid_face = .true.
  endif

!-----------------------------------------------------------------------------------------
  return
  end function Is_on_cuboid_face

!=========================================================================================
  pure function Get_cuboid_face(Point,cell)
!=========================================================================================
! This function tests to which face of the cuboid "cell" the point P (type coordinates)
! belongs to. This function is very similar to function Is_on_cuboid_face but returns
! an array of booleans instead of just one
! 
! Point belongs to the surface of cell is at least one of the condition are met:
! Test if Point is on one of the face.
! If point belongs to cuboid, than the following statements are true:
! d(cell%x(1),cell%x(2)) = d(x,cell%x(1))+d(x,cell%x(2))
! d(cell%y(1),cell%y(2)) = d(y,cell%y(1))+d(y,cell%y(2))
! d(cell%z(1),cell%z(2)) = d(z,cell%z(1))+d(z,cell%z(2))
! Point belongs to one of the face of the cuboid if at least one
! of its coordinates verifies the hyperplane equation.
! Example, If Point belongs to cuboid cell and has y=y1, then we can affirm that
! point belong to the face of equation y=y1 of the cuboid considered.
! Hence we will test that at least one of the d(cell%i,Point(i)) is less than criterion
! Returns false if not
! The function returns an array of logical. The element(s) receiving true are the elements
! of the cuboid for which Point belongs to.
! ----------------------------------------
! The criterion to be on the face is:
! Criterion_face
! Note: Unilike Is_one_cuboid_face, this function does not test if the point is inside
! the cuboid or not. This is done to avoid performance decreases
!----------------------------------------------------------------------------------------

  implicit none

  type(cuboid), intent(in)        :: cell
  type(coordinates), intent(in)   :: Point

  logical, dimension(3,2) :: Get_cuboid_face

  Get_cuboid_face = .false.
! If P is on the surface of the cuboid, it is also inside the cuboid
  
  ! First dimension
  if (metric(Point%x,cell%x(1))<= criterion_face) Get_cuboid_face(1,1) = .true.
  if (metric(Point%x,cell%x(2))<= criterion_face) Get_cuboid_face(1,2) = .true.
! Second dimension
  if (metric(Point%y,cell%y(1))<= criterion_face) Get_cuboid_face(2,1) = .true.
  if (metric(Point%y,cell%y(2))<= criterion_face) Get_cuboid_face(2,2) = .true.
! Third dimension
  if (metric(Point%z,cell%z(1))<= criterion_face) Get_cuboid_face(3,1) = .true.
  if (metric(Point%z,cell%z(2))<= criterion_face) Get_cuboid_face(3,2) = .true.
  
!-----------------------------------------------------------------------------------------
  return
  end function Get_cuboid_face

!=========================================================================================
  pure function Is_inside_cuboid(Point,cell)
!=========================================================================================
! This function tests if the point P (type coordinates) is inside a cuboid "cell" or not
! Returns true if the following conditions are met:
! d(cell%x(1),cell%x(2)) = d(x,cell%x(1))+d(x,cell%x(2))
! d(cell%y(1),cell%y(2)) = d(y,cell%y(1))+d(y,cell%y(2))
! d(cell%z(1),cell%z(2)) = d(z,cell%z(1))+d(z,cell%z(2))
! Use Equality_dim to store quantity: abs(d(cell%x(1),cell%x(2)) - (d(x,cell%x(1))+d(x,cell%x(2))))
! Then do test on all elements of Equality_dim. They should less than epsilon machine
! Returns false if not

  implicit none

  type(cuboid), intent(in)     :: cell
  type(coordinates), intent(in) :: Point

  logical :: Is_inside_cuboid

! Local variables
  real, dimension(3) :: Equality_dim

  Is_inside_cuboid = .true.

  if (Point%x<cell%x(1)) then 
    Is_inside_cuboid = .false.
    return
  endif

  if (Point%x>cell%x(2)) then 
    Is_inside_cuboid = .false.
    return
  endif

  if (Point%y<cell%y(1)) then 
    Is_inside_cuboid = .false.
    return
  endif

  if (Point%y>cell%y(2)) then 
    Is_inside_cuboid = .false.
    return
  endif

  if (Point%z<cell%z(1)) then 
    Is_inside_cuboid = .false.
    return
  endif

  if (Point%z>cell%z(2)) then 
    Is_inside_cuboid = .false.
    return
  endif

!  
!   Equality_dim = 0.0d+0
! 
!   Equality_dim(1) = abs(metric(cell%x(1),cell%x(2))-(metric(Point%x,cell%x(1))+metric(Point%x,cell%x(2))))
! 
!   if (Equality_dim(1).gt.criterion) then 
!     Is_inside_cuboid = .false.
!     return
!   endif
! 
!   Equality_dim(2) = abs(metric(cell%y(1),cell%y(2))-(metric(Point%y,cell%y(1))+metric(Point%y,cell%y(2))))
! 
!   if (Equality_dim(2).gt.criterion) then 
!     Is_inside_cuboid = .false.
!     return
!   endif
! 
!   Equality_dim(3) = abs(metric(cell%z(1),cell%z(2))-(metric(Point%z,cell%z(1))+metric(Point%z,cell%z(2))))
! 
!   if (Equality_dim(3).gt.criterion) then 
!     Is_inside_cuboid = .false.
!     return
!   endif
! 

!   if (maxval(Equality_dim).le.criterion) then
!     Is_inside_cuboid = .true.
!   else
!     Is_inside_cuboid = .false.
!   endif

!  if (Equality_dim(1).gt.criterion) Is_inside_cuboid = .false.
!  if (Equality_dim(2).gt.criterion) Is_inside_cuboid = .false.
!  if (Equality_dim(3).gt.criterion) Is_inside_cuboid = .false.
  

!-----------------------------------------------------------------------------------------
  return
  end function Is_inside_cuboid

!=========================================================================================
  elemental function metric_1d(x1,x2)
!=========================================================================================
! this function returns the distance between the real number x1,x2, which are coordinates
! along a specific dimension of an euclidian space. 
! The distance is given by the relation: metric(x1,x2) = abs(x1-x2)

  implicit none
 
  real, intent(in) :: x1, x2
  real :: metric_1d

  metric_1d = abs(x1-x2) 

!-----------------------------------------------------------------------------------------
  return
  end function metric_1d

!=========================================================================================
  pure function metric_3d(Point_1,Point_2)
!=========================================================================================
! this function returns the Euclidian distance between two point of the 3D space.
! The distance is given by the relation: metric_3d(Point_1,Point_2) = 
! sqrt(Sum_iDim{(Point_2%xi-Point_3%xi)^2})

  implicit none
 
  real, dimension(3), intent(in) :: Point_1, Point_2
  real :: metric_3d

  ! Local Variable

  real, dimension(3) :: buffer
  integer :: i_dim

  metric_3d  = 0.0d+0

  do i_dim = 1, 3
    buffer(i_dim) = metric_1d(Point_2(i_dim),Point_1(i_dim))**2
  enddo

  metric_3d = dsqrt(sum(buffer))

!-----------------------------------------------------------------------------------------
  return
  end function metric_3d

!=========================================================================================
  pure function volume(cell)
!=========================================================================================
! This function returns the volume of a cuboid (rectangular parallelepiped)

  implicit none

  type(cuboid), intent(in) :: cell
 
  real :: volume
  real :: dx1x2, dy1y2, dz1z2 

  volume = 0.0d+0

  dx1x2 = metric(cell%x(1),cell%x(2))
  dy1y2 = metric(cell%y(1),cell%y(2))
  dz1z2 = metric(cell%z(1),cell%z(2))

  volume = dx1x2*dy1y2*dz1z2

!-----------------------------------------------------------------------------------------
  return
  end function volume

!=========================================================================================
  pure function surface_3D(cell)
!=========================================================================================
! This function returns the surface of a cuboid (rectangular parallelepiped)

  implicit none

  type(cuboid), intent(in) :: cell
 
  real :: surface_3D
  real :: dx1x2, dy1y2, dz1z2 

  surface_3D = 0.0d+0

  dx1x2 = metric(cell%x(1),cell%x(2))
  dy1y2 = metric(cell%y(1),cell%y(2))
  dz1z2 = metric(cell%z(1),cell%z(2))

  surface_3D = 2.0d+0*dx1x2*dy1y2 +2.0*dx1x2*dz1z2+2.0*dy1y2*dz1z2

!-----------------------------------------------------------------------------------------
  return
  end function surface_3D

!=========================================================================================
  pure function surface_2D(rectangle)
!=========================================================================================
! This function returns the surface of a facet (rectangular shape)
! A facet has 4 verteces {p1,P2,P3,P4}.
! Its surface is computed by taking the absolute values of the curl of two vectors,
! for example P1P2 and P1P3

  implicit none

  type(facet), intent(in) :: rectangle
 
  real :: surface_2D

  real, dimension(3) :: Vector1, Vector2, Vector3

  surface_2D = 0.0d+0

! Vector1 = defined has vector P1P2
! Vector3 = defined has vector P1P3
! with Pi = rectangle%vertex(i)

  Vector1(1) = rectangle%vertex(2)%x-rectangle%vertex(1)%x
  Vector1(2) = rectangle%vertex(2)%y-rectangle%vertex(1)%y
  Vector1(3) = rectangle%vertex(2)%z-rectangle%vertex(1)%z

  Vector2(1) = rectangle%vertex(3)%x-rectangle%vertex(1)%x
  Vector2(2) = rectangle%vertex(3)%y-rectangle%vertex(1)%y
  Vector2(3) = rectangle%vertex(3)%z-rectangle%vertex(1)%z

  Vector3 = cross_product(Vector1,Vector2)

  surface_2D = dsqrt(DOT_PRODUCT(Vector3,Vector3))

!-----------------------------------------------------------------------------------------
  return
  end function surface_2D

!========================================================================================
  pure function cross_product(Vector1,Vector2)
!========================================================================================

  implicit none

  real, dimension(3), intent(in) :: Vector1, Vector2
  real, dimension(3)             :: cross_product

  cross_product = 0.0d+0

  cross_product(1) = Vector1(2)*Vector2(3)-Vector1(3)*Vector2(2)
  cross_product(2) = Vector1(3)*Vector2(1)-Vector1(1)*Vector2(3)
  cross_product(3) = Vector1(1)*Vector2(2)-Vector1(2)*Vector2(1)

!-----------------------------------------------------------------------------------------
  return 
  end function cross_product

!========================================================================================
  pure function compute_Center(rectangle)
!========================================================================================
! This function computes the the coordinates of a facet barycenter (rectangular shape)
! A facet has 4 verteces {p1,P2,P3,P4}.
! Barycenter is calculated doing:
! compute_Center%x = 0.25*sum(rectangle%vertex(1:4)%x)
! compute_Center%y = 0.25*sum(rectangle%vertex(1:4)%y)
! compute_Center%z = 0.25*sum(rectangle%vertex(1:4)%z)
!----------------------------------------------------------------------------------------

  implicit none

  type(facet), intent(in) :: rectangle
 
  type(coordinates)       :: compute_Center

  compute_Center = Null_coordinates()

  compute_Center%x = 0.25*(sum(rectangle%Vertex(1:4)%x))
  compute_Center%y = 0.25*(sum(rectangle%Vertex(1:4)%y))
  compute_Center%z = 0.25*(sum(rectangle%Vertex(1:4)%z))

!-----------------------------------------------------------------------------------------
  return 
  end function compute_Center

!========================================================================================
  pure function Null_coordinates() result(Null)
!========================================================================================
! This function returns a null structure <coordinates>
!----------------------------------------------------------------------------------------

  implicit none

  type(coordinates) :: Null

  Null%x = 0.0d+0
  Null%y = 0.0d+0
  Null%z = 0.0d+0

!-----------------------------------------------------------------------------------------
  return 
  end function Null_coordinates

!========================================================================================
  function Get_Facet_Coordinates(rectangle) result (Indices)
!========================================================================================
! Returns the coordinates (mx,my,mz) of a facet that is associated with on point of the 
! domain boundaries.
! Argument in is type facet
! Returns Indices = array of integer of dimension 4. The first three are the indices (=[mx, my,mz])
! The last one the index of the plane (1 for South, 2 for North, 3 for Right,
!                                      4 for Left, 5 for Bottom, 6 for Top)

  implicit none

! Variable passed in 
  type(facet), intent(in) :: rectangle

! Variable passed out
  integer, dimension(4)   :: Indices

! Local Variables
  integer :: i_fix_coord

  Indices = 0

! Need to know which boundaries facet belongs to. This is given by the pointer %boundary
! First, check to ensure it is associated with a boundary. If not write error msg.

  if (.NOT.(associated(rectangle%boundary))) call Error_MSG(3)

! Get value of the index corresponding to the boundary equation "i_fix_coord"
! (i.e. if Boundary is South or North then indice = 3,
!       if Boundary is Right or Left then indice = 1, 
!       if Boundary is Bottom or Top then indice = 2, 
! 1 - South : (x,y), eq: z = zz(1)
! 2 - North : (x,y), eq: z = zz(mz)
! 3 - Right : (y,z), eq: x = xx(1)
! 4 - Left  : (y,z), eq: x = xx(mx)
! 5 - Bottom: (x,z), eq: y = yy(1) 
! 6 - Top   : (x,z), eq: y = yy(my)

  i_fix_coord = rectangle%boundary%fixed_dim

! Return the index of the boundary plane.

  Indices(4)  = rectangle%boundary%plane

  Select case (i_fix_coord)
    case(1)
      Indices(1) = rectangle%boundary%coordinate
      Indices(2) = rectangle%i_dim1
      Indices(3) = rectangle%i_dim2
    case(2)
      Indices(1) = rectangle%i_dim1
      Indices(2) = rectangle%boundary%coordinate
      Indices(3) = rectangle%i_dim2
    case(3)
      Indices(1) = rectangle%i_dim1
      Indices(2) = rectangle%i_dim2
      Indices(3) = rectangle%boundary%coordinate
    case default
      call Error_MSG(4)
  end select

!-----------------------------------------------------------------------------------------
  return 
  end function Get_Facet_Coordinates

!=========================================================================================
  pure function Point2Real(Point) result(Point_Real)
!=========================================================================================
! Basic function transforms a type coordinate into a real vector of dimension 3

  implicit none

! Variables passed in
  type(coordinates), intent(in) :: Point

! Variables passed oute
  real, dimension(3) :: Point_Real

! First init Point_real with 0 (precaution)
  Point_Real = 0.0
  
  Point_Real(1) = Point%x
  Point_Real(2) = Point%y
  Point_Real(3) = Point%z

!-----------------------------------------------------------------------------------------
  return 
  end function Point2Real

!=========================================================================================
  pure function Real2Point(Point_Real) result(Point)
!=========================================================================================
! Basic function transforms a real vector of dimension 3 into a type coordinate

  implicit none

! Variables passed in
  real, dimension(3), intent(in) :: Point_Real

! Variables passed out
  type(coordinates):: Point

! First init Point with 0 (precaution)
  Point  = Null_coordinates()
  
  Point%x = Point_Real(1)
  Point%y = Point_Real(2)
  Point%z = Point_Real(3)

!----------------------------------------------------------------------------------------
  return 
  end function Real2Point

!========================================================================================
  pure subroutine Exit_point(Point,cell,direction,Point_exit,i_err)
!========================================================================================
! Point: entry point inside the cuboid <cell>.
! Direction: unitary vector defining the direction of propagation of the ray. Defined 
! with respect to the base assigned with:
! e1 = (Cell%x(2)-Cell%x(1))/metric(Cell%x(2),Cell%x(1))
! e2 = (Cell%y(2)-Cell%y(1))/metric(Cell%y(2),Cell%y(1))
! e3 = (Cell%z(2)-Cell%z(1))/metric(Cell%z(2),Cell%z(1))
!
! If i_err = 0, then subroutine functioned perfectly
! If i_err = 1, then Point is OUTSIDE cell
! If i_err = 2, then direction is null vector
! If i_err = 3, then the solution has not been found
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in

  type(cuboid), intent(in)      :: cell
  type(coordinates), intent(in) :: Point, direction

! Variables passed out

  type(coordinates), intent(out) :: Point_exit
  integer, intent(out)           :: i_err

! Local variables:

  type(coordinates)    :: Point_buffer
  type(coordinates)    :: norm_direction

  real, dimension(3)   :: Point_Real_buffer, Point_Real
  real, dimension(3)   :: norm_Real_direction

  real, dimension(3,2) :: Coordinates_cuboid, d_pp
  real, dimension(3)   :: base_cuboid

  real :: norm, x_prime, y_prime, z_prime, x_p, y_p, z_p, distance_PointPexit

  integer :: i_dim, i_element

  logical, dimension(3,2) :: Is_on_face

  Point_exit = Null_coordinates()

  i_err = 0

! Determine which face of the cuboid the point is on 
! If none, then the problem is irrelevant. Return 
! error message, i_err = 1

  Is_on_face = Get_cuboid_face(Point,cell)
 
  if (.not.(Is_inside_cuboid(Point,cell))) then
    i_err = 1 
    return
  endif

! Check if norm of direction is strictly greater than 0. If not, return i_err = 2
! Point_buffer = origin (0, 0, 0)
  Point_buffer      = Null_coordinates()
  Point_Real_buffer = Point2Real(Point_buffer)
  norm              = metric(Point_Real_buffer,Point2Real(direction))

  if (norm==0.0d+0) then
    i_err = 2
    return
  endif
  
! Normalize direction vector
  if (norm == 1.0d+0) then 
    norm_Real_direction = Point2Real(direction)
  else
    norm_Real_direction = Point2Real(direction)/norm
  endif

! Change type from Real to Point
  norm_direction = Real2Point(norm_Real_direction)

! For each of calculation, Transfer coordinates of the cuboid to the variables coordinates_cuboid.
! These are the coordinates of the vertexe (x1,y1,z1) and (x2,y2,z2)

  Coordinates_cuboid(1,1:2) = cell%x(1:2)
  Coordinates_cuboid(2,1:2) = cell%y(1:2)
  Coordinates_cuboid(3,1:2) = cell%z(1:2)

! local based is oriented along vectors: Coordinates_cuboid(1,2)-Coordinates_cuboid(1,1)
! base_cuboid is introduced to keep consistency with this definition regarless on the value
! of the coordinates

  do i_dim = 1, 3
    base_cuboid(i_dim) = sign(1.0,Coordinates_cuboid(i_dim,2)-Coordinates_cuboid(i_dim,1))
  enddo

! Translates Point into real vector

  Point_real = Point2Real(Point)

! Treat each dimensions sequencially
 
  d_pp = 0.0d+0

!----------------------------------------------------------------------------------------
! Loop over dimension of space
  do i_dim = 1, 3
! Test if the normalized direction vector coordinate along dimension i_dim is null or not
! If null, cycle to other dimension.
    if (metric(norm_Real_direction(i_dim),0.0d+0)<=tiny(norm)) cycle

    do i_element = 1, 2
      if (Is_on_face(i_dim,i_element)) then 
        d_pp(i_dim,i_element) = -1.0d+0
      else
        d_pp(i_dim,i_element) = (Coordinates_cuboid(i_dim,i_element)-Point_real(i_dim)) &
                                /(norm_Real_direction(i_dim)*base_cuboid(i_dim))
      endif
    enddo
  enddo  

! The distance from entry point to the exit point is the smallest strictly positive
  
  distance_PointPexit = MINVAL(d_pp, MASK=d_pp.gt.0.0d+0) 
!  distance_PointPexit = MINVAL(d_pp, MASK=d_pp.gt.criterion)

  Point_exit = Coordinates_translation(Point,Real2Point(norm_Real_direction*base_cuboid),distance_PointPexit)

! ! Final Test: Ensure Point_exit is on the cuboid face. If not, consider the subroutine to
! ! have failed
! 
!   if (.not.Is_on_cuboid_face(Point_exit,cell)) then
!     i_err = 3 
!     Point_exit = Point
!   endif  

!-----------------------------------------------------------------------------------------
  return
  end subroutine Exit_point

!========================================================================================
  subroutine Error_exit_point(i_plane,j_point,i_ray,i_path,ierr,io)
!========================================================================================
! This subroutine interprete the error flag raised by Exit_point and print the
! appropriate statements
!
! Variables in:
! -------------
! i_plane : boundary index
! j_point : point index
! i_ray   : index for ray considered that raised the error flag
! i_path  : index for path segment along the i_ray considered that raised the error flag
! ierr    : error flag raised by Exit_point subroutine
! io      : input output unit
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in 

  integer, intent(in) :: i_plane
  integer, intent(in) :: j_point
  integer, intent(in) :: i_ray
  integer, intent(in) :: i_path
  integer, intent(in) :: ierr
  integer, intent(in) :: io


  if (ierr/=0) then
    write(io,*) 'Error in findpath3d_new!'
    write(io,*) 'For plane: ', i_plane, ', point: ', j_point , '\n'
    write(io,*) 'ray: ', i_ray , ', path: ', i_path , '\n'

    select case(ierr)
      case (1)
        write(io,*) 'Point is OUTSIDE cell.'
      case(2)
        write(io,*) 'Direction is NULL vector.'
      case(3)
        write(io,*) 'Solution has NOT been found.'
      case default
        write(io,*) 'Unknown error.'
    end select
  endif
!-----------------------------------------------------------------------------------------
  return
  end subroutine Error_exit_point

!========================================================================================
  pure function Coordinates_translation(Point,direction,norm) result(Translated_Point)
!========================================================================================
! Compute the new coordinates (Translated_Point) of Point after a translation along
! direction defined by vector direction and of norm nom.
! Coordinates given in the same referentiel

  implicit none

! Variables passed in
  type(coordinates), intent(in) :: Point, direction
  real, intent(in)              :: norm
  
! Variables passed out
  type(coordinates) :: Translated_Point

! Local variables
  real, dimension(3) :: Point_Buffer
  real, dimension(3) :: normalized_direction

  real :: norm_direction

  Translated_Point     = Null_coordinates()
  normalized_direction = 0.0d+0

! Normalize direction vector
  Point_buffer = 0.0d+0
 
  norm_direction = metric(Point_buffer,Point2Real(direction))

  if ((norm_direction/=0.0d+0).and.(norm_direction/=1.0d+0))  then 
      normalized_direction = Point2Real(direction)/norm_direction
  else
      normalized_direction = Point2Real(direction)
  endif
 
  Point_Buffer     = Point2Real(Point) + norm*normalized_direction
  Translated_Point = Real2Point(Point_Buffer)
    
!----------------------------------------------------------------------------------------
  return
  end function Coordinates_translation


!========================================================================================
  function Get_Cell_Index(Point, Index_help) result(Index_Cell)
!========================================================================================
! Returns the indices of the cell containing point "Point"
! Important: This function assumes the cells to have a regular shape. It takes avantages 
! that the cell are rectangular and that there is exactly mx-1 cells in x direction
! my-1 cells in y direction, and mz-1 cells in z direction. 
! If this assumption is not true, then this function needs to be rewritten
! Help can be provided throught the optional argument index_help (integer, dimension(3))
! If provided, the function will test the nearby cells (9 cells maximum)
! If unsuccessful, the function will then proceed with xx, yy, zz
! Of course, this subroutine assumes that xx, yy, zz are sorted in ascending order
! If the indices are all zero, then the Point does not belong to the domain
! This function will write an error message and will stop the execution
!----------------------------------------------------------------------------------------

  implicit none

! Variable pased in
  type(coordinates), intent(in)   :: Point
! Optional variable passed in
  integer, dimension(3), optional :: Index_help
! Variable passed out

  integer, dimension(3) :: Index_Cell

! local variables

  integer :: ix, iy, iz
  integer :: ub_ix, lb_ix
  integer :: ub_iy, lb_iy
  integer :: ub_iz, lb_iz 

  logical :: found_x, found_y, found_z

! Initialize Index_Cell to zero

  Index_Cell = 0

  found_x = .false.
  found_y = .false.
  found_z = .false.

! Test if point is inside the domain (tolerance is machine epsilon)
! If not print error message and leave

  if (.not.(Is_inside_cuboid(Point,Global_Cell))) call Error_MSG(5)

  if (present(Index_help)) then 
    ! Define local upper and lower bound for local search near Index_help
    ub_ix = mx - 1
    ub_iy = my - 1
    ub_iz = mz - 1

    lb_ix = 1
    lb_iy = 1
    lb_iz = 1

    ! Make sure Index_help elements are between the bounds:
    ! 2 < Index_help(1)<=mx-1
    ! 1 < Index_help(2)<=my-1
    ! 1<=Index_help(3)<=mz-1
    if (1 < Index_help(1)  )  lb_ix = Index_help(1)-1
    if (Index_help(1)< mx-1)  ub_ix = Index_help(1)+1

    if (1 < Index_help(2)  )  lb_iy = Index_help(2)-1
    if (Index_help(2)< my-1)  ub_iy = Index_help(2)+1

    if (1 < Index_help(3)  )  lb_iz = Index_help(3)-1
    if (Index_help(3)< mz-1)  ub_iz = Index_help(3)+1
  
    ! Loop in x 
 
    do ix = lb_ix, ub_ix
      if ((xx(ix)<=Point%x).and.(Point%x<=xx(ix+1))) then
        Index_Cell(1) = ix
        found_x = .true.
        exit
      else
        cycle
      endif
    enddo    

    do iy = lb_iy, ub_iy
      if ((yy(iy)<=Point%y).and.(Point%y<=yy(iy+1))) then
        Index_Cell(2) = iy
        found_y = .true.
        exit
      else
        cycle
      endif
    enddo    

    do iz = lb_iz, ub_iz
      if ((zz(iz)<=Point%z).and.(Point%z<=zz(iz+1))) then
        Index_Cell(3) = iz
        found_z = .true.
        exit
      else
        cycle
      endif
    enddo        

  endif ! (PRESENT(Index_help)) conditions

!----------------------------------------------------------------------------------------
! Get indices of Point. If found_<x,y,z> then skip appropriate part

  if (.not.(found_x)) then
! Loop over x direction. Find ix such as xx(ix)<=P%x<xx(ix+1)
    do ix = 1, mx-1
      if ((xx(ix)<=Point%x).and.(Point%x<=xx(ix+1))) then
        Index_Cell(1) = ix
        found_x = .true.
        exit
      else
        cycle
      endif
    enddo
  endif


  if (.not.(found_y)) then
! Loop over y direction. Find iy such as yy(iy)<=P%y<yy(iy+1)
    do iy = 1, my-1
      if ((yy(iy)<=Point%y).and.(Point%y<=yy(iy+1))) then
        Index_Cell(2) = iy
        found_y = .true.
        exit
      else
        cycle
      endif
    enddo
  endif

  if (.not.(found_z)) then
! Loop over z direction. Find iz such as zz(iz)<=P%z<zz(iz+1)
    do iz = 1, mz-1
      if ((zz(iz)<=Point%z).and.(Point%z<=zz(iz+1))) then
        Index_Cell(3) = iz
        found_z = .true.
        exit
      else
        cycle
      endif
    enddo
  endif

!----------------------------------------------------------------------------------------
! Verification; Test that Point is in cell_tot(ix,iy,iz) with
! ix = Index(1), iy = Index(2), iz = Index(2)

  ix = Index_Cell(1)
  iy = Index_Cell(2)
  iz = Index_Cell(3)

  if (.not.(Is_inside_cuboid(Point,cell_tot(ix,iy,iz)))) then
    write(*,*) 'Problem in Get_Cell_Index'
    Index_Cell = 0
  endif

!----------------------------------------------------------------------------------------
  return
  end function Get_Cell_Index

!========================================================================================
  function Get_Facet_Index(Point) result(Indices)
!========================================================================================
! Returns the indices of the facets containing point "Point"
!
! Of course, this subroutine assumes that xx, yy, zz are sorted in ascending order
! If the first indice is zero, then the Point does not belong to the domain.
! Note: This function returns the indeces of the closest facets in case Point does not
! belong to the boundary. This was done to account for round-off errors that can become
! problematic. 
! This function will write an error message and will stop the execution
! Indices has three element:
! Indices(1) : Index of the boundary containing the facet: 1 South, 2 North, 3 Right
! 4 Left, 5 Bottom, 6 Top
! Note: The tolerance is higher in this function to account for point that are slightly
! outside the domain.
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in:
  type(coordinates), intent(in) :: Point

! Variables passed out:
  real, dimension(3) :: Indices 

! Local variables
  real, dimension(:), pointer :: Dim1, Dim2

  real, dimension(3) :: Point_Real 

  integer            :: ix, iy, iz, i_boundary, i_dim, i_dim1, i_dim2
  integer            :: n_dim1, n_dim2, index_dim1, index_dim2

! Initialize Index_Cell to zero
  Indices = 0

! Test if point is on the surface of the cuboid that the domain defines
! If not print error message and leave
! DEPRECATED: if (.not.(Is_on_cuboid_face(Point,Global_Cell))) call Error_MSG(5)
! DO NOT TEST IN THE POINT BELONG TO THE CUBOID FACE. THE REASON IS THAT THE TOLERANCE
! IN Is_on_cuboid_face IS VERY LOW, AS IT SHOULD BE. HOWEVER THE ORIGIN POINT OF THE RAY
! SHOULD HAVE HIGHER TOLERANCE

! Loop over Boundaries and test if metric(Fixed_coordinate,Point) is less than criterion

  Point_Real = Point2Real(Point)

! Find first indices
  do i_boundary = 1, 6

    i_dim = Boundaries(i_boundary)%fixed_dim
    select case(i_boundary)
      case (1,3,5)
        if (Point_real(i_dim)<=Boundaries(i_boundary)%Fixed_coordinate+criterion_face) then
          Indices(1) = i_boundary
          exit
        endif
      case (2,4,6)
        if (Point_real(i_dim)>=Boundaries(i_boundary)%Fixed_coordinate-criterion_face) then
          Indices(1) = i_boundary
          exit
        endif
    end select
  enddo

  i_boundary = Indices(1)

! Incomplete. Stopped here on 06/06/13. Need to check if code is okay

  n_dim1 = size(Boundaries(i_boundary)%element,1)
  n_dim2 = size(Boundaries(i_boundary)%element,2)

  if (associated(Dim1)) nullify(Dim1)
  if (associated(Dim2)) nullify(Dim2)

  Select case (i_boundary)
    case(1:2)
      index_dim1 = 1
      index_dim2 = 2
      Dim1 => xx
      Dim2 => yy
    case(3:4)
      index_dim1 = 2
      index_dim2 = 3
      Dim1 => yy
      Dim2 => zz
    case(5:6)
      index_dim1 = 1
      index_dim2 = 3
      Dim1 => xx
      Dim2 => zz
    case default
      call Error_MSG(4)
  end select

!------------------------------------------------------------------------------
! Find indices for first dimension
  if (Point_Real(index_dim1)<Dim1(1))    then
    Indices(2) = 1
  elseif(Point_Real(index_dim1)>Dim1(n_dim1)) then
    Indices(2) = n_dim1
  else
    do i_dim1 = 1, n_dim1-1
      if ((Dim1(i_dim1)<=Point_Real(index_dim1)).and.(Point_Real(index_dim1)<=Dim1(i_dim1+1))) then
        Indices(2) = i_dim1
        exit
      else
        cycle
      endif
    enddo
  endif

!------------------------------------------------------------------------------
! Find indices for second dimension
  if (Point_Real(index_dim2)<Dim2(1))    then
    Indices(3) = 1
  elseif(Point_Real(index_dim2)>Dim2(n_dim2)) then
    Indices(3) = n_dim2
  else
    do i_dim2 = 1, n_dim2-1
      if ((Dim2(i_dim2)<=Point_Real(index_dim2)).and.(Point_Real(index_dim2)<=Dim2(i_dim2+1))) then
        Indices(3) = i_dim2
        exit
      else
        cycle
      endif
    enddo
  endif

!----------------------------------------------------------------------------------------
  return
  end function Get_Facet_Index

!=========================================================================================
  subroutine findpath3d_new(i_plane,io)
!=========================================================================================
! This subroutine tracks the path of each ray that fired from each facets of the boundary
! designated by the integer plane.
! This (private) subroutine update the public variable ray(:,:) (type beam).
! Recall: ray(n_rays,n_points). Initialized in subroutine <allocate_radiation_arrays>
! 
! Arguments in:
!--------------
! integer :: i_plane ! index of the boundary to treat: 
!                  ! Recall : 
!                   1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top 
! io      :: unit of the output file
!
! Arguments out:
!---------------
! NONE
!
! Recall:
!--------
!   type beam
!     type(facet), pointer                :: Facet_origin
!     type(cuboid), dimension(:), pointer :: Passed_CV   ! get adresses of control volume crossed
!     real,  allocatable, dimension(:)    :: Length      ! beam length of ray in the control volume
!     integer :: n_path    
!     integer :: iray2
!     integer :: iphi
!     integer :: itheta
!     integer :: edge_final
!   end type beam
!
!-----------------------------------------------------------------------------------------

  implicit none  

  integer, intent(in)    :: i_plane
  integer, intent(inout) :: io

! local declarations
 
  type(beam)            :: ray_buffer
  type(coordinates)     :: P0, P1, P01 

  type(coordinates) :: vec_dir1  ! coordinates of the direction vector 
                                 ! expressed in based attached to boundaries(i_plane)

  type(coordinates) :: vec_dir   ! coordinates of the direction vector in "absolute" base
                                 ! i.e. base associated with South plane: (x,y,z) 

  real, parameter   :: criterion_3 = 2.5*criterion ! small increment by 2.5 time machine precision
                                                   ! needed to assess which next cell is crossed
                                                   ! recall that criterion for belonging to cell or
                                                   ! facet is criterion

  real :: cos_theta, sin_theta        ! Save value of cosine and sine of theta  
  real :: cos_phi, sin_phi            ! Save value of cosine and sine of phi
  real :: lambda, lambda_buffer       ! Coefficient for point increment

  integer :: l, m, i_path

  integer, pointer :: i, j, k

  integer, dimension(3) :: Index_cell, Index_facet, I_help
  integer, target :: i_outer, i_inner
  integer :: n_outer, n_inner

  integer :: i_phi, i_theta
  integer :: i_ray, j_point              ! Dummy indices for ray tables
  integer :: n_max                       ! Greatest number of cell crossed: = sqrt([mx-1]^2+[my-1]^2+[mz-1]^2)
  integer :: ierr


  logical :: end_ray

! end of declarations
!----------------------------------------------------------------------------------------  

! Estimate the maximum number of cells crossed:
! Length of longuest path is given by sqrt(dim_x^2+dim_y^2+dim_z^2). Assume dim_i = Ni,i={x,y,x} 
! Add extra points for long rays. Problems can arise otherwise
  n_max = ceiling(dsqrt(float((mx+3)**2+(my+3)**2+(mz+3)**2))) + max(mx,my,mz)
  
!----------------------------------------------------------------------------------------
! Loop over each point of the boundary i_plane
! Recall:
!
!------------------------------------
! Plane   | outer loop | inner loop |
! 1-South |  x         |     y      |
! 2-North |  x         |     y      |
! 3-Right |  y         |     z      |
! 4-Left  |  y         |     z      |
! 5-Bottom|  z         |     x      |
! 6-Top   |  z         |     x      |
!------------------------------------
!
! Important considerations for the determination of j_point

  select case(i_plane)
    case(1:2)
      n_outer = mx-1
      n_inner = my-1

      i => i_outer
      j => i_inner
      k => Boundaries(i_plane)%coordinate
    case(3:4)
      n_outer = my-1
      n_inner = mz-1

      i => Boundaries(i_plane)%coordinate
      j => i_outer
      k => i_inner
    case(5:6)
      n_outer = mz-1
      n_inner = mx-1

      i => i_inner
      j => Boundaries(i_plane)%coordinate
      k => i_outer
  end select

!----------------------------------------------------------------------------------------
! Loop over all the facets of boundaries(i_plane)

OUTERLOOP: do i_outer = 1, n_outer  ! loop over outer loop points

INNERLOOP: do i_inner = 1, n_inner  ! loop over inner loop points

      j_point = Get_Points_Index(i,j,k,i_plane)

!----------------------------------------------------------------------------------------
! Loop over rays.    
! Initialized i_ray

      i_ray = 0

LOOPTHETA: do i_theta = 1, n_theta    ! loop on polar angle theta

        cos_theta = dcos(theta(i_theta))
        sin_theta = dsin(theta(i_theta))

LOOPPHI: do i_phi = 1, n_phi(i_theta) ! loop on azimutal angle

          i_ray = i_ray + 1  

          cos_phi = dcos(phi(i_theta)%azimuth(i_phi))
          sin_phi = dsin(phi(i_theta)%azimuth(i_phi))

          ! Create direction vector associated with local base
          ! vec_dir1

          vec_dir1%x = cos_phi*sin_theta
          vec_dir1%y = sin_phi*sin_theta
          vec_dir1%z = cos_theta

          ! Compute coordinates of vec_dir1 in South basis

          vec_dir%x=DOT_PRODUCT(Change_plane(i_plane,1,1,1:3),Point2Real(vec_dir1))
          vec_dir%y=DOT_PRODUCT(Change_plane(i_plane,1,2,1:3),Point2Real(vec_dir1))
          vec_dir%z=DOT_PRODUCT(Change_plane(i_plane,1,3,1:3),Point2Real(vec_dir1))

          ! Create ray_buffer and initialize it.
          ! Eventually ray(i_ray,i_point) = ray_buffer

          if (allocated(ray_buffer%Passed_CV)) deallocate(ray_buffer%Passed_CV)
          if (allocated(ray_buffer%Length))    deallocate(ray_buffer%Length)

! TO REMOVE AFTER DEBUGGING
          if (allocated(ray_buffer%Point))    deallocate(ray_buffer%Point)
! END
          allocate (ray_buffer%Passed_CV(n_max))
          allocate (ray_buffer%Length(n_max))   ; ray_buffer%Length = 0

! TO REMOVE AFTER DEBUGGING
          allocate (ray_buffer%Point(n_max+1))
! END

          ray_buffer%n_path     = 1
          ray_buffer%iray2      = 0
    
          if (associated(ray_buffer%Facet_origin)) nullify(ray_buffer%Facet_origin)

          ray_buffer%iphi       = i_phi
          ray_buffer%itheta     = i_theta
          ray_buffer%edge_final = i_plane

          ! Initialize counter for cells crossed
          i_path = 1

          ! Warning. Boundaries(i_plane)%element(1:mx-1,1:mz-1)
          ! (and NOT Boundaries(i_plane)%element(1:mz-1,1:mx-1)
          ! Get coordinates of the center of the facet. Store it in P0
 
          if ((i_plane == 5).or.(i_plane == 6)) then
            ray_buffer%Passed_CV(i_path)%cub => Boundaries(i_plane)%element(i_inner,i_outer)%cell
            P0                               =  Boundaries(i_plane)%element(i_inner,i_outer)%center
          else
            ray_buffer%Passed_CV(i_path)%cub => Boundaries(i_plane)%element(i_outer,i_inner)%cell 
            P0                               =  Boundaries(i_plane)%element(i_outer,i_inner)%center
          endif

          P01     = P0
          end_ray = .false.

!       TO REMOVE EVENTUALLY
          ray_buffer%Point(i_path) = P0

!----------------------------------------------------------------------------------------
! Loop until end of ray is reached, i.e. origin of the ray has been found
ALONGRAY: do while((.not.end_ray).and.(i_path<=n_max))

            ! Compute exit point P1. Use P01 instead P0 as a precaution. Save Length.
            ! Test if any boundary has been reached
            call Exit_point(P01,ray_buffer%Passed_CV(i_path)%cub,vec_dir,P1,ierr)
           
            if (i/=0) call Error_exit_point(i_plane,j_point,i_ray,i_path,ierr,io)

            ! Save length P0P1 (length of the ray that intersects the cell
            ray_buffer%Length(i_path) = metric(Point2Real(P0),Point2Real(P1))
            
            ! Update P0 using value of P1
            P0 = P1
            ray_buffer%Point(i_path+1) = P1

! ! Debugging code
!             if ((i_ray == 2).and.(j_point==87)) then 
!               print *, 'Here is P0: ', P0%x, P0%y, P0%z
!               print *, 'Here is P1: ', P1%x, P1%y, P1%z
!               print *, ' '
!               print *, 'Here is Global_Cell: x1 = ', Global_Cell%x(1), ' x2 = ', Global_Cell%x(2)
!               print *, 'Here is Global_Cell: y1 = ', Global_Cell%y(1), ' y2 = ', Global_Cell%y(2)
!               print *, 'Here is Global_Cell: z1 = ', Global_Cell%z(1), ' z2 = ', Global_Cell%z(2)
!               print *, ' '
!               print *, 'Results Is_on_cuboid_face(P1,Global_Cell): ', Is_on_cuboid_face(P1,Global_Cell)
!               print *, 'Results Is_inside_cuboid(P1,Global_Cell) : ', Is_inside_cuboid(P1,Global_Cell)
!               print *, 'Nmax: ', n_max, ' , here is i_path ', i_path
!               print *, '=========================================='
!             endif
! ! End debugging code

            ! Test if P1 is on a boundary face. Use variable Global_Cell

            if ((.not.(Is_on_cuboid_face(P1,Global_Cell))).and.(Is_inside_cuboid(P1,Global_Cell))) then 

              call RANDOM_NUMBER(lambda)  ! Coefficient for increment
              lambda = lambda*1.0d+4

              ! Calculate P01 which is given by P1+vec_dir*3*criterion
              P01 = Coordinates_translation(P1,vec_dir,lambda*criterion_3)

              ! Give some help based on previous knowledge of ray_buffer%Passed_CV(i_path)
              I_help(1)  = ray_buffer%Passed_CV(i_path)%cub%i_x
              I_help(2)  = ray_buffer%Passed_CV(i_path)%cub%i_y
              I_help(3)  = ray_buffer%Passed_CV(i_path)%cub%i_z

              ! Find the cell that contains P01              
              Index_cell = Get_Cell_Index(P01,I_help)
              
              !--------------------------------------------------------------------------
              ! Test to ensure that PO1 is not on the cuboid face 
              ! (Problem with associated cell can arise otherwise)
              ! If P01 is on the cuboid face, then multiply the increment by factor lambda
              ! compute new P01 and test again. Lambda follows an exponential progression
              ! Needed. Otherwise bad rays can be constructed

              do while (Is_on_cuboid_face(P01,cell_tot(Index_cell(1),Index_cell(2),Index_cell(3))) &
                        .and.(lambda<3.0d+4)) 

                call RANDOM_NUMBER(lambda_buffer) ! Variable stride to improve out of cuboid_face condition

                lambda     = lambda + lambda_buffer*1.0d+4
                P01        = Coordinates_translation(P1,vec_dir,lambda*criterion_3)
                Index_cell = Get_Cell_Index(P01,I_help) ! Find the cell that contains P01. Give help to reduce cost
                I_help     = Index_cell
              enddo

              i_path     = i_path + 1
              ray_buffer%Passed_CV(i_path)%cub => cell_tot( Index_cell(1),Index_cell(2),Index_cell(3) )

            else
              ! The exit point is on a boundary.
              ! Get facet that contains P1
              end_ray                 = .true.
              Index_facet             = Get_Facet_Index(P1)
              ray_buffer%n_path       = i_path
              ray_buffer%Facet_origin => Boundaries(Index_facet(1))%element(Index_facet(2),Index_facet(3))

            endif

          enddo ALONGRAY ! (while condition, exit when reach end of ray)

!----------------------------------------------------------------------------------------
          ! Save the ray into global variable
          ! Allocate memory for ray
          if (allocated(ray(i_ray,j_point)%Passed_CV)) deallocate(ray(i_ray,j_point)%Passed_CV )
          if (allocated(ray(i_ray,j_point)%Length))    deallocate(ray(i_ray,j_point)%Length    )

! TO REMOVE EVENTUALLY. FOR DEBUGGING PURPOSES ONLY
          if (allocated(ray(i_ray,j_point)%Point))     deallocate(ray(i_ray,j_point)%Point     )
          allocate (ray(i_ray,j_point)%Point(ray_buffer%n_path + 1))
! END

          allocate (ray(i_ray,j_point)%Passed_CV(ray_buffer%n_path))
          allocate (ray(i_ray,j_point)%Length(ray_buffer%n_path))

          ! Assign values of ray_buffer into ray
          ray(i_ray,j_point)%Passed_CV    =  ray_buffer%Passed_CV(1:ray_buffer%n_path)
          ray(i_ray,j_point)%Length       =  ray_buffer%Length(1:ray_buffer%n_path)
          ray(i_ray,j_point)%Facet_origin => ray_buffer%Facet_origin

          ray(i_ray,j_point)%Point        =  ray_buffer%Point(1:(ray_buffer%n_path+1))

          ray(i_ray,j_point)%n_path       =  ray_buffer%n_path
          ray(i_ray,j_point)%iray2        =  0 ! This is assigned later

          ray(i_ray,j_point)%iphi         =  ray_buffer%iphi
          ray(i_ray,j_point)%itheta       =  ray_buffer%itheta
          ray(i_ray,j_point)%edge_final   =  ray_buffer%edge_final

        enddo  LOOPPHI ! End loop azimutal angle
      enddo LOOPTHETA  ! End loop polar angle theta
  
    enddo INNERLOOP    ! End inner loop
  enddo  OUTERLOOP     ! End outer loop

!----------------------------------------------------------------------------------------
  return
  end subroutine findpath3d_new


!========================================================================================
  subroutine Update_Cell(yspecies,pressure,temp,io)
!========================================================================================
! This function updates the following attribute of the global variable cell_tot
! cell_tot%c_abs       : Planck mean absorption coefficient of the cell
! cell_tot%Temp        : Temperature of the cellule
! cell_tot%steradiancy : BlackBody radiant flux emitted. Units in W/m^2W/str.
!----------------------------------------------------------------------------------------
  implicit none

! Variables passed in
  real, intent(in), dimension(:,:,:)   :: temp      ! Temperature, UNITS = Kelvin
  real, intent(in), dimension(:,:,:)   :: pressure  ! Pressure,    UNITS = Pascal
  real, intent(in), dimension(:,:,:,:) :: yspecies  ! Species mass fraction. UNITS = N/A

  integer, intent(in), optional :: io               ! i/o unit

! Local variables

  integer :: ix, iy, iz

! Update absorption first 
  call get_absorption(yspecies,pressure,temp) 

! Update cell Temperature second (important: do it before computing the blackbody steradiancy
  cell_tot%Temp = 0.0
  call Transf_Ptwise_CV(temp,cell_tot%Temp)

! Update black body radiant lux emitted

  do iz = 1, mz-1
    do iy = 1, my-1
      do ix = 1, mx-1
        cell_tot(ix,iy,iz)%steradiancy = sigma/pi*cell_tot(ix,iy,iz)%Temp**4
      enddo
    enddo
  enddo

!----------------------------------------------------------------------------------------
  return
  end subroutine Update_Cell

!========================================================================================
  subroutine get_absorption(yspecies,pressure,temp,io)
!========================================================================================
! First compute the absorption coefficients for each grid point (kappa(nx,ny,nz))
! and then average it to get that for the CV                    (cell_tot%c_abs(mx-1,my-1,mz-1))
! An alternative way is to first get the temperature, Ys etc for the CV 
! and then compute absorption coefficient based on them
! However, this requires more data transfer
!----------------------------------------------------------------------------------------

!  use chemistry_m, only : molwt_c

  implicit none

! Variables passed in
  real, intent(in), dimension(:,:,:)   :: temp      ! Temperature, UNITS = Kelvin
  real, intent(in), dimension(:,:,:)   :: pressure  ! Pressure,    UNITS = Pascal
  real, intent(in), dimension(:,:,:,:) :: yspecies  ! Species mass fraction. UNITS = N/A

  integer, intent(in), optional :: io               ! i/o unit

! local declarations
! kappa_soot=C_soot*fv*T [1/m] , note: Yoo use 1117.97*fv*T
  real, parameter :: C_soot = 1298.0

  real, allocatable, dimension(:,:,:,:) :: x_mol
  real, allocatable, dimension(:,:,:,:) :: dkp
!  real sum

  integer, dimension(3) :: dim_xyz
  integer :: n_species, nx, ny ,nz

!  integer :: i, j, k, isc, jj, kk
!----------------------------------------------------------------------------------------
! Ensure the variable have same dimensions

  if (all(shape(temp).NE.shape(pressure))) then
    call Error_MSG(1)

  elseif (all(shape(temp).NE.shape(yspecies(:,:,:,1)))) then
    call Error_MSG(2)
  else
    dim_xyz = shape(temp)
    nx = dim_xyz(1)
    ny = dim_xyz(2)
    nz = dim_xyz(3)
    n_species = size(yspecies,DIM=4)
  endif

!----------------------------------------------------------------------------------------
! Allocate some space

  allocate(x_mol(nx,ny,nz,n_species))
  allocate(dkp(nx,ny,nz,n_species))

  x_mol = 0.0
  dkp   = 0.0
!----------------------------------------------------------------------------------------
! Compute kappa at grid points

!  where (temp >= 2.0d+3)
!    kappa = 1.0d+6
!  elsewhere
!    kappa = 5.0d-2
!  end where
   kappa = 1.0d-2
!   do k = 1, nz
!     do j = 1, ny
!       do i = 1, nx
! 
!         sum = 0.0
!         kappa(i,j,k) = 0.0
!  CVL TEST5
!        kappa(i,j,k) = 0.01 
! Compute species mole fraction

!         do isc = 1, nsc+1
!           x_mol(isc) = yspecies(i,j,k,isc) * molwt_c(isc)
!           sum = sum + x_mol(isc)
!         enddo
! 
!         do isc = 1, nsc+1
!           x_mol(isc) = x_mol(isc) / sum
!         enddo
! !------------------------------------------------------------------------------
! ! Compute coefficients
! ! Recall: i_count_rad defined in opthinrad_m module, represents number of 
! ! species radiation data set in opthinrad.in file
! 
!         do jj = 1, i_count_rad
! 
! ! Compute Planck mean absorption coefficients [1/(atm*m)]
! ! TNF data: 
! ! dkp(j)= as(j,1) + as(j,2)*(1000/T) + as(j,3)*(1000/T)^2 + as(j,4)*(1000/T)^3 
! !        + as(j,5)*(1000/T)^4 + as(j,6)*(1000/T)^5
! 
!           dkp(jj) = as(jj,6)
! 
!           do kk = 5, 1, -1
!             dkp(jj) = dkp(jj)*(1.0d+3/temp(i,j,k)) + as(jj,kk)
!           enddo
! 
! ! Compute total mean absorption coefficient [1/m]
! 
!           if(index_rad(jj).gt.0) then
!             kappa(i,j,k) = kappa(i,j,k) +                                      & 
!                          x_mol(index_rad(jj))*pressure(i,j,k)/pres_atm*dkp(jj)
!           endif
! 
!         enddo
!     
! ! Add the effect of soot
!         if (i_soot==1) then
!           kappa(i,j,k) = kappa(i,j,k) + C_soot*fv(i,j,k)*temp(i,j,k)
!         endif 

!       enddo
!     enddo
!   enddo !

!---------------------------------------------------------------------------------------
! get c_abs from kappa. Call subroutine Transf_Ptwise
!------------------------------------------------------------------------------------------

  call Transf_Ptwise_CV(kappa, cell_tot%c_abs)

  deallocate(x_mol)
  deallocate(dkp)
!---------------------------------------------------------------------------------------
  return
  end subroutine get_absorption

!=========================================================================================
  subroutine init_rad_geometry(nx,ny,nz,x,y,z,io)
!=========================================================================================
! prepare beam lengh and path for each ray fired from boundary
!-----------------------------------------------------------------------------------------
  use topology_m, only : myid

  implicit none

! declaration passed in
  integer :: io
  integer :: nx,ny,nz
  real :: x(nx),y(ny),z(nz)

!-----------------------------------------------------------------------------------------
! track and record the path for every ray
! order boundary points as : south (x,y), north (x,-y), right (y,z), left (-y,z),
!                            bottom (z,x), top (z,-x)
! plane normal: south (z), north (-z), right (x), left (-x), bottom (y), top (-y)
! order rays fired : counter-clockwise for every boundary point 
! ------------------ It is so for tracking rays from ghost cells

! South
  call findpath3d_new(1,io)

! North
  call findpath3d_new(2,io)
  
! Right plane: (y,z,x)
  call findpath3d_new(3,io)

! Left plane: (-y,z,-x)
  call findpath3d_new(4,io)

! Bottom plane (z,x,y)
  call findpath3d_new(5,io)

! Top plane (z,-x,-y)
  call findpath3d_new(6,io)

  call Find_Associated_Ray(ray)

! CVL: Print specific ray  call print_ray(1, 1806)
! write header
  if(myid.eq.0) then
    write(io,*) 'Finish ray geometry construction.'
  endif    

  return
  end subroutine init_rad_geometry
!----------------------------------------------------------------------------------------

!=========================================================================================
  subroutine initialize_ChangePlane
!=========================================================================================
! Subroutine checked by Vivien Lecoustre on 07/05/12
! Status: Valid
! Transformation Matrices P_ij - USED IN DTM SUBROUTINE
! (south) -> (north)  (x,y,z) -> (x,-y,-z)
!       |1  0  0|   
! P_SN =|0 -1  0|  
!       |0  0 -1| 

! (south) -> (right)  (x,y,z) -> (y,z,x)
!       |0  1  0|   
! P_SR =|0  0  1|  
!       |1  0  0| 

! (south) -> (left) (x,y,z) -> (y,-z,-x)
!       |0 -1  0|   
! P_SL =|0  0  1|  
!       |-1 0  0| 

! (south) -> (bottom) (x,y,z) -> (z,x,y)
!       |0  0  1|   
! P_SB =|1  0  0|  
!       |0  1  0| 

! (south) -> (top)  (x,y,z) -> (z,-x,-y)
!       |0  0  1|   
! P_ST =|-1 0  0|  
!       |0 -1  0| 
!-----------------------------------From North Plane-------------------------------------
! (north) -> (south)  (x,-y,-z) -> (x,y,z)
!       |1  0  0|   
! P_NS =|0 -1  0| (=P_SN) 
!       |0  0 -1| 

! (north) -> (right)  (x,-y,-z) -> (y,z,x)
!       |0 -1  0|   
! P_NR =|0  0 -1|  
!       |1  0  0| 

! (north) -> (left) (x,-y,-z) -> (-y,z,-x)
!       |0  1  0|   
! P_NL =|0  0 -1|   
!       |-1 0  0| 

! (north) -> (bottom) (x,-y,-z) -> (z,x,y)
!       |0  0 -1|   
! P_NB =|1  0  0|  
!       |0 -1  0| 

! (north) -> (top)  (x,-y,-z) -> (z,-x,-y)
!       |0  0 -1|   
! P_NT =|-1 0  0|  
!       |0  1  0| 
!--------------------------From Right Plane----------------------------------------------
! (right) -> (south)  (y,z,x) -> (x,y,z)
!       |0  0  1|   
! P_RS =|1  0  0| (=P_SB)  
!       |0  1  0| 

! (right) -> (north)  (y,z,x) -> (x,-y,-z)
!       |0  0  1|   
! P_RN =|-1 0  0| (=P_ST) 
!       |0 -1  0| 

! (right) -> (left) (y,z,x) -> (-y,z,-x)
!       |-1 0  0|   
! P_RL =|0  1  0|   
!       |0  0 -1| 

! (right) -> (bottom) (y,z,x) -> (z,x,y)
!       |0  1  0|   
! P_RB =|0  0  1| (=P_SR)
!       |1  0  0| 

! (right) -> (top)  (y,z,x) -> (z,-x,-y)
!       |0  1  0|   
! P_RT =|0  0 -1| (=P_NL)  
!       |-1 0  0| 

!--------------------------From Left Plane----------------------------------------------
! (left) -> (south)  (-y,z,-x) -> (x,y,z)
!       |0  0 -1|   
! P_LS =|-1 0  0|  (=P_NT)
!       |0  1  0| 

! (left) -> (north)  (-y,z,-x) -> (x,-y,-z)
!       |0  0 -1|   
! P_LN =|1  0  0| (=P_NB)
!       |0 -1  0| 

! (left) -> (right) (-y,z,-x) -> (y,z,x)
!       |-1 0  0|   
! P_LR =|0  1  0|  (=P_RL)
!       |0  0 -1| 

! (left) -> (bottom) (-y,z,-x) -> (z,x,y)
!       |0  1  0|   
! P_LB =|0  0 -1| (=P_NL)
!       |-1 0  0| 

! (left) -> (top)  (-y,z,-x) -> (z,-x,-y)
!       |0  1  0|   
! P_LT =|0  0  1| (=P_SR)  
!       |1  0  0| 

!--------------------------From Bottom Plane----------------------------------------------
! (bottom) -> (south)  (z,x,y) -> (x,y,z)
!       |0  1  0|   
! P_BS =|0  0  1|  (=P_SR)   
!       |1  0  0| 

! (bottom) -> (north)  (z,x,y) -> (x,-y,-z)
!       |0  1  0|   
! P_BN =|0  0 -1|  (=P_NL)
!       |-1 0  0| 

! (bottom) -> (right) (z,x,y) -> (y,z,x)
!       |0  0  1|   
! P_BR =|1  0  0|  (=P_SB)
!       |0  1  0| 

! (bottom) -> (left) (z,x,y) -> (-y,z,-x)
!       |0  0 -1|   
! P_BL =|1  0  0| (=P_NB)
!       |0 -1  0| 

! (bottom) -> (top)  (z,x,y) -> (z,-x,-y)
!       |1  0  0|   
! P_LT =|0 -1  0| (=P_SN)  
!       |0  0 -1| 

!--------------------------From Top Plane----------------------------------------------
! (top) -> (south)  (z,-x,-y) -> (x,y,z)
!       |0 -1  0|   
! P_TS =|0  0 -1|  (=P_NR)   
!       |1  0  0| 

! (top) -> (north)  (z,-x,-y) -> (x,-y,-z)
!       |0 -1  0|   
! P_TN =|0  0  1| (=P_SL)
!       |-1 0  0| 

! (top) -> (right) (z,-x,-y) -> (y,z,x)
!       |0  0 -1|   
! P_TR =|1  0  0|  (=P_NB)
!       |0 -1  0| 

! (top) -> (left) (z,-x,-y) -> (-y,z,-x)
!       |0  0  1|   
! P_TL =|1  0  0|  (=P_SB) 
!       |0  1  0| 

! (top) -> (bottom)  (z,-x,-y) -> (z,x,y)
!       |1  0  0|   
! P_TB =|0 -1  0| (=P_SN)  
!       |0  0 -1| 
!
!  Recall: 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top
!----------------------------------------------------------------------------------------
  implicit none

! Variable transformation matrices. Will populate Change_plane(6,6,3,3)
  real, dimension(3,3) :: P_SS, P_SN, P_SR, P_SL, P_SB, P_ST
  real, dimension(3,3) :: P_NS, P_NN, P_NR, P_NL, P_NB, P_NT
  real, dimension(3,3) :: P_RS, P_RN, P_RR, P_RL, P_RB, P_RT
  real, dimension(3,3) :: P_LS, P_LN, P_LR, P_LL, P_LB, P_LT
  real, dimension(3,3) :: P_BS, P_BN, P_BR, P_BL, P_BB, P_BT
  real, dimension(3,3) :: P_TS, P_TN, P_TR, P_TL, P_TB, P_TT

  Change_plane = 0.0

  P_SS = 0.0
  P_SN = 0.0
  P_SR = 0.0
  P_SL = 0.0
  P_SB = 0.0
  P_ST = 0.0

  P_NN = 0.0
  P_NS = 0.0
  P_NR = 0.0
  P_NL = 0.0
  P_NB = 0.0
  P_NT = 0.0

  P_RR = 0.0
  P_RS = 0.0
  P_RN = 0.0
  P_RL = 0.0
  P_RB = 0.0
  P_RT = 0.0

  P_LL = 0.0
  P_LS = 0.0 
  P_LN = 0.0
  P_LR = 0.0 
  P_LB = 0.0 
  P_LT = 0.0

  P_BB = 0.0
  P_BS = 0.0 
  P_BN = 0.0 
  P_BR = 0.0 
  P_BL = 0.0 
  P_BT = 0.0

  P_TT = 0.0
  P_TS = 0.0 
  P_TN = 0.0
  P_TR = 0.0 
  P_SL = 0.0 
  P_TB = 0.0

! Assigne values to matrices
  P_SS(1,1) = 1.0
  P_SS(2,2) = 1.0
  P_SS(3,3) = 1.0

  P_SN(1,1) = 1.0
  P_SN(2,2) = -1.0
  P_SN(3,3) = -1.0

  P_SR(1,2) = 1.0
  P_SR(2,3) = 1.0
  P_SR(3,1) = 1.0

  P_SL(1,2) = -1.0
  P_SL(2,3) = 1.0
  P_SL(3,1) = -1.0

  P_SB(1,3) = 1.0
  P_SB(2,1) = 1.0
  P_SB(3,2) = 1.0

  P_ST(1,3) = 1.0
  P_ST(2,1) = -1.0
  P_ST(3,2) = -1.0

  P_NN = P_SS
  P_NS = P_SN

  P_NR(1,2) = -1.0
  P_NR(2,3) = -1.0
  P_NR(3,1) = 1.0

  P_NL(1,2) = 1.0
  P_NL(2,3) = -1.0
  P_NL(3,1) = -1.0

  P_NB(1,3) = -1.0
  P_NB(2,1) = 1.0
  P_NB(3,2) = -1.0

  P_NT(1,3) = -1.0
  P_NT(2,1) = -1.0
  P_NT(3,2) =  1.0

  P_RR = P_SS
  P_RS = P_SB
  P_RN = P_ST

  P_RL(1,1) = -1.0
  P_RL(2,2) =  1.0
  P_RL(3,3) = -1.0

  P_RB = P_SR
  P_RT = P_NL

  P_LS = P_NT 
  P_LN = P_NB
  P_LR = P_RL 
  P_LB = P_NL 
  P_LT = P_SR
  P_LL = P_SS
  
  P_BS = P_SR 
  P_BN = P_NL
  P_BR = P_SB
  P_BL = P_NB
  P_BT = P_SN
  P_BB = P_SS

  P_TS = P_NR 
  P_TN = P_SL
  P_TR = P_NB
  P_TL = P_SB
  P_TB = P_SN
  P_TT = P_SS

  Change_plane(1,1,:,:) = P_SS
  Change_plane(1,2,:,:) = P_SN
  Change_plane(1,3,:,:) = P_SR
  Change_plane(1,4,:,:) = P_SL
  Change_plane(1,5,:,:) = P_SB
  Change_plane(1,6,:,:) = P_ST

  Change_plane(2,1,:,:) = P_NS
  Change_plane(2,2,:,:) = P_NN
  Change_plane(2,3,:,:) = P_NR
  Change_plane(2,4,:,:) = P_NL
  Change_plane(2,5,:,:) = P_NB
  Change_plane(2,6,:,:) = P_NT

  Change_plane(3,1,:,:) = P_RS
  Change_plane(3,2,:,:) = P_RN
  Change_plane(3,3,:,:) = P_RR
  Change_plane(3,4,:,:) = P_RL
  Change_plane(3,5,:,:) = P_RB
  Change_plane(3,6,:,:) = P_RT

  Change_plane(4,1,:,:) = P_LS
  Change_plane(4,2,:,:) = P_LN
  Change_plane(4,3,:,:) = P_LR
  Change_plane(4,4,:,:) = P_LL
  Change_plane(4,5,:,:) = P_LB
  Change_plane(4,6,:,:) = P_LT

  Change_plane(5,1,:,:) = P_BS
  Change_plane(5,2,:,:) = P_BN
  Change_plane(5,3,:,:) = P_BR
  Change_plane(5,4,:,:) = P_BL
  Change_plane(5,5,:,:) = P_BB
  Change_plane(5,6,:,:) = P_BT

  Change_plane(6,1,:,:) = P_TS
  Change_plane(6,2,:,:) = P_TN
  Change_plane(6,3,:,:) = P_TR
  Change_plane(6,4,:,:) = P_TL
  Change_plane(6,5,:,:) = P_TB
  Change_plane(6,6,:,:) = P_TT

!----------------------------------------------------------------------------------------
  return
  end subroutine initialize_ChangePlane

!=========================================================================================
  subroutine Transf_Ptwise_CV(Ptwise_var, Cv_var)
!=========================================================================================
! Create a Control Volume variable based on the informations provided by a pointwise 
! variables.
! INPUT
!------
!
! Ptwise_var : pointwise variable to be transformed in Control Volume variables.
!              Type: Real. Dimensions(nx,zny,nz)
!
! OUTPUT
!-------
!
! Cv_var : Control volume variable. 
!	   Type: Real. Dimensions (mx,my,mz)
!
! This module has been developped to handle 3D cases
!----------------------------------------------------------------------------------------
! 
!                    ______________________
!                   /|                    /|
!                  / |  neighbor(6)      / |
!                 /  |  Plan(North)     /  |
!                /   |                 /   |
!               /_____________________/    |
!              |     |                |    |
!              |     |          neighb|r(4)|
!              |     |            plan|top)|
! neighbor(1)  |     |                |    | neighbor(2) 
! Plan(right)  |     |________________|____| Plan(left)
!              |    /                 |    /
!              |   /                  |   /
!           neighbor(3)               |  /
!           Plan(bottom)              | /
!              |/_____________________|/
!
!                   neighbor(5)
!                   Plan(south)
!
!         (z)^  (y)
!            | /
!            |/--> (x)
! 
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in  
  real, dimension(:,:,:), intent(in) :: Ptwise_var
  
! Variables passed out
  real, dimension(:,:,:), intent(out) :: Cv_var

! Local Variables
  real, allocatable, dimension(:,:,:) :: Ptwise_var_ghost 

  integer, dimension(3) :: dim_xyz

  integer :: ncvx, ncvy, ncvz  ! size of Cv_var in each direction

  integer i, j, k

! Get Cv_var size in each direction
  dim_xyz = shape(Cv_var)

  ncvx = dim_xyz(1)
  ncvy = dim_xyz(2)
  ncvz = dim_xyz(3)

! Initialization

  Cv_var = 0.0

!----------------------------------------------------------------------------------------
! Get Variable Ptwise_var appended with its needed ghost planes and edges in new Variable
! Ptwise_var_ghost

! Deallocate Ptwise_var_ghost. Its dimensions will be set automatically by subroutine
! Exch_ghost_plan

  if(allocated(Ptwise_var_ghost)) deallocate(Ptwise_var_ghost)

  call Exch_ghost_plan(Ptwise_var, Ptwise_var_ghost)

!----------------------------------------------------------------------------------------
! Create Control Volume variables

  do k = 1, ncvz
    do j = 1, ncvy
      do i = 1, ncvx
        Cv_var(i,j,k)=0.125*(Ptwise_var_ghost(i,j  ,k  ) + Ptwise_var_ghost(i+1,j  ,k  )+&
                             Ptwise_var_ghost(i,j+1,k  ) + Ptwise_var_ghost(i+1,j+1,k  )+&
                             Ptwise_var_ghost(i,  j,k+1) + Ptwise_var_ghost(i+1,j  ,k+1)+&
                             Ptwise_var_ghost(i,j+1,k+1) + Ptwise_var_ghost(i+1,j+1,k+1))
      enddo
    enddo
  enddo
!----------------------------------------------------------------------------------------
! Free space

  if(allocated(Ptwise_var_ghost)) deallocate(Ptwise_var_ghost) 

  return
!----------------------------------------------------------------------------------------
  end subroutine Transf_Ptwise_CV

!========================================================================================
  subroutine Exch_ghost_plan(Variable,Variable_with_Ghost)
!========================================================================================
! Variables. Assumed shape
! INPUT
!------
! Variable : type real, array (:,:,:). Either pointwise or cell. Dimensions Known
! 
! OUTPUT
!-------
! Variable_with_Ghost : type real, array (:,:,:) Either pointwise or cell. Corresponds
!                                                to Variable added with ghost planes and 
!                                                edges if necessary.
!----------------------------------------------------------------------------------------
  include 'mpif.h'

! Variables passed in  
  real, dimension(:,:,:), intent(in) :: Variable
  
! Variables passed out
  real, allocatable, dimension(:,:,:), intent(out) :: Variable_with_Ghost

! Local Variables
  integer, dimension(3) :: dim_xyz

  integer :: nx,  ny,  nz         ! size of Variable in each direction
  integer :: nx1, ny1, nz1        ! size of Variable_with_Ghost in each direction

! Coordinates of the process that has the memory windows
  
  integer, dimension(3) :: coordinates
  integer               :: target_id ! Rank of remote process

! Define new MPI_DATA_TYPE
  integer :: XY_plane, YZ_plane, XZ_plane
  integer :: XY_plane_ghost, YZ_plane_ghost, XZ_plane_ghost

  integer :: Z_vector, Y_vector, X_vector              
  integer :: Z_vector_ghost, Y_vector_ghost, X_vector_ghost
              
! ! Planes for ghost cells
!   real, allocatable, dimension(:,:) :: Var_left
!   real, allocatable, dimension(:,:) :: Var_north
!   real, allocatable, dimension(:,:) :: Var_top
!  
! ! Edges for ghost cells
!   real, allocatable, dimension(:) :: Var_ed_north
!   real, allocatable, dimension(:) :: Var_ed_top
!   real, allocatable, dimension(:) :: Var_ed_left
!  
! ! Vertex in x(nx+1),y(ny+1),z(nz+1) 
!   real :: Var_corner

! Windows objects
  integer :: win_var                ! Windows object for Variable

  integer(KIND=MPI_ADDRESS_KIND) :: sizereal, target_disp

  sizereal = size_type_real
! Get Variable size in each direction

  dim_xyz = shape(Variable)

  nx = dim_xyz(1)
  ny = dim_xyz(2)
  nz = dim_xyz(3)

! Create the Momery Window and return the window oject

  call MPI_Win_Create(Variable,nx*ny*nz*sizereal,size_type_real,MPI_INFO_NULL,cart_comm,win_var,ierr_MPI_rad)

! Note 3D. Need seven patchs of information
! 3 planes: [x(nx+1)], [y(ny+1)], [z(nz+1)]
! 3 edges: [x(nx+1),y(ny+1)], [x(nx+1),z(ny+1)], [y(nx+1),z(ny+1)]
! 1 corner point [x(nx+1), y(ny+1), z(nz+1)]

! !----------------------------------------------------------------------------------------
! ! Allocate variables needed for path of informations
! ! Planes
!   if (allocated(Var_left))  deallocate(Var_left)
!   if (allocated(Var_top))   deallocate(Var_top)
!   if (allocated(Var_north)) deallocate(Var_north)
! 
!   allocate(Var_left(ny,nz))
!   allocate(Var_top(nx,nz))
!   allocate(Var_north(nx,ny))
! 
! ! Edges
! 
!   if (allocated(Var_ed_left))  deallocate(Var_ed_left)
!   if (allocated(Var_ed_north)) deallocate(Var_ed_north)
!   if (allocated(Var_ed_top))   deallocate(Var_ed_top)
!  
!   allocate(Var_ed_left(ny))
!   allocate(Var_ed_north(nx))
!   allocate(Var_ed_top(nz))
! 
!   Var_left     = 0.0
!   Var_north    = 0.0
!   Var_top      = 0.0
!  
!   Var_ed_north = 0.0
!   Var_ed_top   = 0.0
!   Var_ed_left  = 0.0
!  
!   Var_corner   = 0.0

! Define dimension of Variable_with_Ghost (nx1,ny1,nz1)
! If ghost_<i> is true for direction <i>=(x,y,z) then n<i>1 = n<i>+1, else n<i>1 = n<i>

  if (ghost_x) then 
    nx1 = nx+1
  else
    nx1 = nx
  endif

  if (ghost_y) then 
    ny1 = ny+1
  else
    ny1 = ny
  endif

  if (ghost_z) then 
    nz1 = nz+1
  else
    nz1 = nz
  endif

! Allocate memory for Variable_with_Ghost

  if(allocated(Variable_with_Ghost)) deallocate(Variable_with_Ghost)
  allocate(Variable_with_Ghost(nx1,ny1,nz1))

! Define MPI_TYPE_VECTOR
! These MPI_TYPE are defined to read Variable of Dimensions(nx,ny,nz)
! YZ_plane, XZ_plane, XY_plane
! X_vector, Y_vector, Z_vector

  call MPI_Type_vector(   nz,     1, nx*ny, TypeReal, Z_vector, ierr_MPI_rad)
  call MPI_Type_vector(   ny,     1,    nx, TypeReal, Y_vector, ierr_MPI_rad)
  call MPI_Type_Contiguous( nx,             TypeReal, X_vector, ierr_MPI_rad)

  call MPI_Type_commit(Z_vector,ierr_MPI_rad)
  call MPI_Type_commit(Y_vector,ierr_MPI_rad)
  call MPI_Type_commit(X_vector,ierr_MPI_rad)

  call MPI_Type_vector(ny*nz,     1,    nx, TypeReal, YZ_plane, ierr_MPI_rad)
  call MPI_Type_vector(   nz,    nx, nx*ny, TypeReal, XZ_plane, ierr_MPI_rad)
  call MPI_Type_contiguous(nx*ny,           TypeReal, XY_plane, ierr_MPI_rad)

  call MPI_Type_commit(YZ_plane,ierr_MPI_rad)
  call MPI_Type_commit(XZ_plane,ierr_MPI_rad)
  call MPI_Type_commit(XY_plane,ierr_MPI_rad)

! These MPI_TYPE are defined to read Variable_with_Ghost of Dimensions(nx1,ny1,nz1)
! YZ_plan_ghost,   XZ_plan_ghost,  XY_plan_ghost
! X_vector_ghost, Y_vector_ghost, Z_vector_ghost

  call MPI_Type_vector(nz, 1, nx1*ny1, TypeReal, Z_vector_ghost, ierr_MPI_rad)
  call MPI_Type_vector(ny, 1,     nx1, TypeReal, Y_vector_ghost, ierr_MPI_rad)
  call MPI_Type_contiguous(nx,         TypeReal, X_vector_ghost, ierr_MPI_rad)

  call MPI_Type_commit(Z_vector_ghost,ierr_MPI_rad)
  call MPI_Type_commit(Y_vector_ghost,ierr_MPI_rad)
  call MPI_Type_commit(X_vector_ghost,ierr_MPI_rad)

!  call MPI_Type_vector(ny*nz, 1,      nx1, sizereal, YZ_plan_ghost, ierr_MPI_rad)

  call MPI_Type_create_hvector(nz, 1,nx1*ny1*sizereal, Y_vector_ghost, YZ_plane_ghost, ierr_MPI_rad) ! Heterogeneous stride
  call MPI_Type_vector(nz, nx, nx1*ny1, TypeReal,       XZ_plane_ghost, ierr_MPI_rad)
  call MPI_Type_vector(ny, nx,     nx1, TypeReal,       XY_plane_ghost, ierr_MPI_rad)

  call MPI_Type_commit(YZ_plane_ghost,ierr_MPI_rad)
  call MPI_Type_commit(XZ_plane_ghost,ierr_MPI_rad)
  call MPI_Type_commit(XY_plane_ghost,ierr_MPI_rad)

! Transfer Variable inner points into Variable_with_Ghost 

  Variable_with_Ghost(1:nx,1:ny,1:nz) = Variable

! Start Remote memory access Get operation
  call MPI_WIN_FENCE(0,win_var,ierr_MPI_rad)

!----------------------------------------------------------------------------------------
! Ghost plane: x-direction 

  if (ghost_x) then
! Get Var_left

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1)+1, myid_cart(2), myid_cart(3)]
! Get its address using MPI_CART_RANK
    coordinates      = myid_cart
    coordinates(1)   = myid_cart(1) +1
    target_disp      = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_left,ny*nz,sizereal,target_id,sizereal,1,YZ_plane,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(nx1,1,1),1,YZ_plane_ghost,target_id,target_disp,1,YZ_plane,win_var,ierr_MPI_rad)

  endif

! Ghost plane: y-direction 

  if (ghost_y) then
! Get Var_top

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1), myid_cart(2)+1, myid_cart(3)]
! Get its address using MPI_CART_RANK
    coordinates      = myid_cart
    coordinates(2)   = myid_cart(2)+1
    target_disp      = 0
    
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_top,nx*nz,sizereal,target_id,sizereal,1,XZ_plane,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(1,ny1,1),1,XZ_plane_ghost,target_id,target_disp,1,XZ_plane,win_var,ierr_MPI_rad)

  endif

! Ghost plane: z-direction 

  if (ghost_z) then
! Get Var_north

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1), myid_cart(2), myid_cart(3)+1]
! Get its address using MPI_CART_RANK
    coordinates      = myid_cart
    coordinates(3)   = myid_cart(3)+1
    target_disp      = 0

    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_north,nx*ny,sizereal,target_id,sizereal,1,XY_plan,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(1,1,nz1),1,XY_plane_ghost,target_id,target_disp,1,XY_plane,win_var,ierr_MPI_rad)

  endif

!----------------------------------------------------------------------------------------
! Ghost edges: Vector along Z-direction

  if (ghost_x.and.ghost_y) then
! Get Var_ed_top

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1)+1, myid_cart(2)+1, myid_cart(3)]
! Get its address using MPI_CART_RANK
    coordinates(1:2) = myid_cart(1:2)+1
    coordinates(3)   = myid_cart(3)
    target_disp      = 0
    
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_ed_top,nz,sizereal,target_id,sizereal,1,Z_vector,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(nx1,ny1,1),1,Z_vector_ghost,target_id,target_disp,1,Z_vector,win_var,ierr_MPI_rad)

  endif

! Ghost edges: Vector along Y-direction 

  if (ghost_x.and.ghost_z) then
! Get Var_ed_left

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1)+1, myid_cart(2), myid_cart(3)+1]
! Get its address using MPI_CART_RANK
    coordinates(1)   = myid_cart(1)+1
    coordinates(2)   = myid_cart(2)
    coordinates(3)   = myid_cart(3)+1
    target_disp      = 0
    
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_ed_left,ny,sizereal,target_id,sizereal,1,Y_vector,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(nx1,1,nz1),1,Y_vector_ghost,target_id,target_disp,1,Y_vector,win_var,ierr_MPI_rad)

  endif

! Ghost edges: Vector along X-direction 

  if (ghost_y.and.ghost_z) then
! Get Var_ed_north

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1), myid_cart(2)+1, myid_cart(3)+1]
! Get its address using MPI_CART_RANK
    coordinates(1)   = myid_cart(1)
    coordinates(2)   = myid_cart(2)+1
    coordinates(3)   = myid_cart(3)+1
    target_disp      = 0
    
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_ed_north,nx,sizereal,target_id,sizereal,1,X_vector,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(1,ny1,nz1),1,X_vector_ghost,target_id,target_disp,1,X_vector,win_var,ierr_MPI_rad)

  endif

!----------------------------------------------------------------------------------------
! Vertex in nx1,ny1,nz1. 
  if (ghost_x.and.ghost_y.and.ghost_z) then
! Get Var_corner

! Get the address of the process where memory will be read from
! In Cartesian, its address is : [myid_cart(1)+1, myid_cart(2)+1, myid_cart(3)+1]
! Get its address using MPI_CART_RANK
    coordinates(1)   = myid_cart(1)+1
    coordinates(2)   = myid_cart(2)+1
    coordinates(3)   = myid_cart(3)+1
    target_disp      = 0
    
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)

!    call MPI_Get(Var_corner,1,sizereal,target_id,0,1,sizereal,win_var,ierr_MPI_rad)
    call MPI_Get(Variable_with_Ghost(nx1,ny1,nz1),1,TypeReal,target_id,target_disp,1,TypeReal,win_var,ierr_MPI_rad)

  endif

! Stop Remote memory access Get operation
  call MPI_WIN_FENCE(0,win_var,ierr_MPI_rad)

! Free memory
  call MPI_Win_Free(win_var,ierr_MPI_rad)

! Free MPI DATATYPE
  call MPI_Type_Free (YZ_plane, ierr_MPI_rad)
  call MPI_Type_Free (XY_plane, ierr_MPI_rad)
  call MPI_Type_Free (XZ_plane, ierr_MPI_rad)

  call MPI_Type_Free (YZ_plane_ghost, ierr_MPI_rad)
  call MPI_Type_Free (XY_plane_ghost, ierr_MPI_rad)
  call MPI_Type_Free (XZ_plane_ghost, ierr_MPI_rad)

  call MPI_Type_Free (X_vector, ierr_MPI_rad)
  call MPI_Type_Free (Y_vector, ierr_MPI_rad)
  call MPI_Type_Free (Z_vector, ierr_MPI_rad)

  call MPI_Type_Free (X_vector_ghost, ierr_MPI_rad)
  call MPI_Type_Free (Y_vector_ghost, ierr_MPI_rad)
  call MPI_Type_Free (Z_vector_ghost, ierr_MPI_rad)

  return
!----------------------------------------------------------------------------------------
  end subroutine Exch_ghost_plan

!========================================================================================
  subroutine Radiative_emission(temp,incident_flux,bc_emission,optional_parameter)
!========================================================================================
! Variables. Assumed shape
! INPUT
!------
! temp          : type real, array (:,:,:). Temperature, units are in Kelvin
! incident_flux : type real, array(:).    Incident radiative heat flux from boundary point
! optional_parameter : type real, optional. parameter for user implementation (ex: time)
! 
! OUTPUT
!-------
! bc_emission : type real,vector. Radiative flux emitted from a boundary point. Units W/m2
!----------------------------------------------------------------------------------------
! Compute the radiative heat emited in one direction from the boundary points.
! Assume Lambertian surface.
! Recall that the points are located at the interface between the Control Volumes and the 
! boundaries
! Need to get temperatures from the ghost planes
! This function requires the knowledge of mx, my, and mz
!----------------------------------------------------------------------------------------

  implicit none

! Variables that are passed in and out
  real, dimension(:,:,:), intent(in)  :: temp
  real, dimension(:), intent(in)      :: incident_flux
  real, dimension(:), intent(out)     :: bc_emission

  real, optional :: optional_parameter ! optional parameter for user implementation

! local variables

! Temperature at a given boundary points
  real :: temp_point

! Temperature array with ghost planes
  real, allocatable, dimension(:,:,:) :: temp_with_ghost

  integer :: i_points

  integer :: i, j, k
  integer :: i_origin

  integer :: nx1, ny1, nz1
  integer :: nx2, nx2b

  if (allocated(temp_with_ghost)) deallocate(temp_with_ghost)
  call Exch_ghost_plan(temp,temp_with_ghost)

! Get dimensions for temp, incident_flux, and emission
! Assumes array lower bounds start at 1.

  nx1 = UBOUND(temp_with_ghost,1)
  ny1 = UBOUND(temp_with_ghost,2)
  nz1 = UBOUND(temp_with_ghost,3)

  nx2  = UBOUND(bc_emission,1)
  nx2b = UBOUND(incident_flux,1)

  if (nx2/=nx2b) print *, 'Problem, dimensions bc_emission and incident_flux should be identical'

  call Update_facet_Temperature(temp_with_ghost) ! Refresh facet temperature

  if(allocated(temp_with_ghost)) deallocate(temp_with_ghost)

  bc_emission = 0.0

!----------------------------------------------------------------------------------------
! get the emissive power of each boundary cell (wall or virtual boundary). The distinction
! will be done at the RTE integration level

XY_PLANES: do i_origin = 1, 2

    if (i_origin==1) then        ! South plane
      k = 1
    elseif (i_origin == 2) then  ! North plane
      k = nz1
    endif

   ! Loop over elements of the Boundary
    do i = 1, nx1-1
      do j = 1, ny1-1

        i_points = Get_Points_Index(i,j,k,i_origin)

        ! Black body emission
        bc_emission(i_points) = Boundaries(i_origin)%emissivity*sigma* & 
                                Boundaries(i_origin)%element(i,j)%Temp**4/pi 
      ! Contribution of the reflected incident radiative flux
        bc_emission(i_points) = bc_emission(i_points) + &
                           (1.0-Boundaries(i_origin)%emissivity)*incident_flux(i_points)/pi 

      enddo  
    enddo
  enddo XY_PLANES 
 
YZ_PLANES: do i_origin = 3, 4

    if (i_origin==3) then        ! Right plane
      i = 1
    else                         ! Left plane
      i = nx1
    endif

   ! Loop over elements of the Boundary
    do j = 1, ny1-1
      do k = 1, nz1-1
      
        i_points = Get_Points_Index(i,j,k,i_origin)
      
      ! Black body emission
        bc_emission(i_points) = Boundaries(i_origin)%emissivity*sigma* &
                              Boundaries(i_origin)%element(j,k)%Temp**4/pi 
      ! Contribution of the reflected incident radiative flux
        bc_emission(i_points) = bc_emission(i_points) + & 
                              (1.0-Boundaries(i_origin)%emissivity)*incident_flux(i_points)/pi

      enddo
    enddo
  enddo YZ_PLANES

XZ_PLANES: do i_origin = 5, 6

   if (i_origin==5) then        ! Bottom plane
     j = 1
   else                         !  Top plane
     j = ny1
   endif

   ! Loop over elements of the Boundary
    do k = 1, nz1-1
      do i = 1, nx1-1
   
        i_points = Get_Points_Index(i,j,k,i_origin)

        ! Black body emission
        bc_emission(i_points) = Boundaries(i_origin)%emissivity*sigma* &
                              Boundaries(i_origin)%element(i,k)%Temp**4/pi 
        ! Contribution of the reflected incident radiative flux
        bc_emission(i_points) = bc_emission(i_points) + &
                              (1.0-Boundaries(i_origin)%emissivity)*incident_flux(i_points)/pi

      enddo
    enddo
  enddo XZ_PLANES

!----------------------------------------------------------------------------------------
  return
  end subroutine Radiative_emission

!========================================================================================
  subroutine Update_facet_Temperature(temp_with_ghost,optional_parameter)
!======================================================================================== 
! This function computes the temperature at the center of each facet that 
! comprises the domain
! It was extracted from the subroutine Radiative_emission and put as its own for increased
! control over the temperature. This is important because it controls the temperature
! at boundary. This function should only be calles within Radiative_emission.
!  Recall: 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-TopPoint
!----------------------------------------------------------------------------------------

  implicit none

  ! Variables passed in

  real, dimension(:,:,:), intent(in) :: temp_with_ghost
  real, optional :: optional_parameter ! optional parameter for user implementation

  ! Local variables

  integer :: i, j, k
  integer :: i_origin

  integer :: nx1, ny1, nz1

  real, dimension(6) :: BC_Temp

  ! Note: If the user wants to implement non-fixed BC for radiation, do it here  

  ! Implement BC_Temp
  BC_Temp(1) =  Temp_South
  BC_Temp(2) =  Temp_North
  BC_Temp(3) =  Temp_Right
  BC_Temp(4) =  Temp_Left
  BC_Temp(5) =  Temp_Bottom
  BC_Temp(6) =  Temp_Top

  ! Get dimension of temp_with_ghost

  nx1 = UBOUND(temp_with_ghost,1)
  ny1 = UBOUND(temp_with_ghost,2)
  nz1 = UBOUND(temp_with_ghost,3)

  ! Loop over planes

  ! South and North planes first
XY_PLANES: do i_origin = 1, 2

    ! If Virtual boundary, use gas temperature (important since it will vary)
    ! If NOT virtual boundary, use value in BC_Temp and cycle
    ! Recall Virtuall boundary is set by the MPI subroutine

    if (.not.is_vb(i_origin)) then
      Boundaries(i_origin)%element(:,:)%Temp = BC_Temp(i_origin)
      cycle
    endif

    if (i_origin==1) then        ! South plane
      k = 1
    elseif (i_origin == 2) then  ! North plane
      k = nz1
    endif

   ! Loop over elements of the Boundary
    do i = 1, nx1-1
      do j = 1, ny1-1

        Boundaries(i_origin)%element(i,j)%Temp = 0.25*( temp_with_ghost(i  ,j  ,k )  + &
                                                        temp_with_ghost(i+1,j  ,k )  + &
                                                        temp_with_ghost(i  ,j+1,k )  + &
                                                        temp_with_ghost(i+1,j+1,k ) )

      enddo  
    enddo
  enddo XY_PLANES 

  ! Right and Left planes first
YZ_PLANES: do i_origin = 3, 4

    ! If Virtual boundary, use gas temperature (important since it will vary)
    ! If NOT virtual boundary, use value in BC_Temp and cycle
    ! Recall Virtuall boundary is set by the MPI subroutine

    if (.not.is_vb(i_origin)) then
      Boundaries(i_origin)%element(:,:)%Temp = BC_Temp(i_origin)
      cycle
    endif

    if (i_origin==3) then        ! Right plane
      i = 1
    else                         ! Left plane
      i = nx1
    endif

   ! Loop over elements of the Boundary
    do j = 1, ny1-1
      do k = 1, nz1-1

        Boundaries(i_origin)%element(j,k)%Temp = 0.25*( temp_with_ghost(i ,j  ,k  )  + &
                                                        temp_with_ghost(i ,j+1,k  )  + & 
                                                        temp_with_ghost(i ,j  ,k+1)  + &
                                                        temp_with_ghost(i ,j+1,k+1) )

      enddo
    enddo
  enddo YZ_PLANES

  ! Bottom and Top
XZ_PLANES: do i_origin = 5, 6

    ! If Virtual boundary, use gas temperature (important since it will vary)
    ! If NOT virtual boundary, use value in BC_Temp and cycle
    ! Recall Virtuall boundary is set by the MPI subroutine

    if (.not.is_vb(i_origin)) then
      Boundaries(i_origin)%element(:,:)%Temp = BC_Temp(i_origin)
      cycle
    endif

   if (i_origin==5) then        ! Bottom plane
     j = 1
   else                         !  Top plane
     j = ny1
   endif

   ! Loop over elements of the Boundary
    do k = 1, nz1-1
      do i = 1, nx1-1
   
        Boundaries(i_origin)%element(i,k)%Temp = 0.25*( temp_with_ghost(i  ,j ,k  )  + & 
                                                        temp_with_ghost(i  ,j ,k+1)  + &
                                                        temp_with_ghost(i+1,j ,k  )  + &
                                                        temp_with_ghost(i+1,j ,k+1) )

      enddo
    enddo
  enddo XZ_PLANES

!----------------------------------------------------------------------------------------
  return
  end subroutine Update_facet_Temperature


!========================================================================================
  subroutine exchange_RIVB
!========================================================================================
! This routine exchanges the Ray Intensity tables at Virtual Boundaries
! Use of Cartesian coordinates accounts implicitly for the treatment of periodic_bc
!----------------------------------------------------------------------------------------
  implicit none

  include 'mpif.h'
  
! Local variables for MPI
  integer :: RIB_XY, RIB_YZ, RIB_XZ
  
  integer,dimension(3) :: coordinates
  integer :: target_id 

! Windows objects
  integer :: win_north, win_south  ! Windows object RIVB along Z direction
  integer :: win_left,  win_right  ! Windows object RIVB along X direction
  integer :: win_top,   win_bottom ! Windows object RIVB along Y direction

  integer(KIND=MPI_ADDRESS_KIND) :: sizereal
  integer(KIND=MPI_ADDRESS_KIND) :: n_points
  integer(KIND=MPI_ADDRESS_KIND) :: target_disp

  target_disp = 0
  sizereal    = size_type_real

! Create MPI derived datatype  
  call MPI_TYPE_CONTIGUOUS(n_rays*(mx-1)*(my-1),TypeReal,RIB_XY,ierr_MPI_rad)
  call MPI_TYPE_CONTIGUOUS(n_rays*(my-1)*(mz-1),TypeReal,RIB_YZ,ierr_MPI_rad)
  call MPI_TYPE_CONTIGUOUS(n_rays*(mx-1)*(mz-1),TypeReal,RIB_XZ,ierr_MPI_rad)

  call MPI_Type_commit(RIB_XY,ierr_MPI_rad)
  call MPI_Type_commit(RIB_YZ,ierr_MPI_rad)
  call MPI_Type_commit(RIB_XZ,ierr_MPI_rad)

! Define remote memory window and create window objects
! Memory window for Left and Right exchange
  n_points = n_rays*(my-1)*(mz-1)*sizereal
  call MPI_Win_Create(rivb_r   ,n_points,size_type_real,MPI_INFO_NULL,cart_comm,win_right,ierr_MPI_rad) ! For MPI_put
  call MPI_Win_Create(rib%Right,n_points,size_type_real,MPI_INFO_NULL,cart_comm,win_left ,ierr_MPI_rad) ! For MPI_get

! Memory window for North and South exchange
  n_points = n_rays*(mx-1)*(my-1)*sizereal
  call MPI_Win_Create(rivb_s   ,n_points,size_type_real,MPI_INFO_NULL,cart_comm,win_south,ierr_MPI_rad) ! For MPI_put
  call MPI_Win_Create(rib%South,n_points,size_type_real,MPI_INFO_NULL,cart_comm,win_north,ierr_MPI_rad) ! For MPI_get

! Memory window for Top and Bottom exchange
  n_points = n_rays*(mx-1)*(mz-1)*sizereal
  call MPI_Win_Create(rivb_b    ,n_points,size_type_real,MPI_INFO_NULL,cart_comm,win_bottom,ierr_MPI_rad) ! For MPI_put
  call MPI_Win_Create(rib%Bottom,n_points,size_type_real,MPI_INFO_NULL,cart_comm,win_top   ,ierr_MPI_rad) ! For MPI_get

! Remember that virtual cartesian topology implicitly accounts 
! for periodic boundary conditions. Use Logical ghost_x, ghost_y, ghost_z to determine 
! if tranfer via MPI_put and MPI_Get is needed.

! Activate communication

  call MPI_WIN_FENCE(0,win_right ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_left  ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_south ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_north ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_bottom,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_top   ,ierr_MPI_rad)

! Transfer between Left and Right
!--------------------------------
  if ((ghost_x).or.(periodic_bc(1))) then
    ! Get the rank of the process of coordinates [X+1,Y,Z]
    coordinates      = myid_cart
    coordinates(1)   = myid_cart(1) + 1
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)
    call MPI_Put(RIB%Left(1,1,1), 1, RIB_YZ, target_id,target_disp,1,RIB_YZ,win_right,ierr_MPI_rad)
    call MPI_Get(rivb_l,          1, RIB_YZ, target_id,target_disp,1,RIB_YZ,win_left ,ierr_MPI_rad)
  endif

! Transfer between Top and Bottom
!--------------------------------
  if ((ghost_y).or.(periodic_bc(2))) then 
    ! Get the rank of the process of coordinates [X,Y+1,Z]
    coordinates      = myid_cart
    coordinates(2)   = myid_cart(2) + 1
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)
    call MPI_Put(RIB%Top(1,1,1), 1, RIB_XZ, target_id,target_disp,1,RIB_XZ,win_bottom,ierr_MPI_rad)
    call MPI_Get(rivb_t,         1, RIB_XZ, target_id,target_disp,1,RIB_XZ,win_top ,ierr_MPI_rad)
  endif

! Transfer between North and South
!---------------------------------
  if ((ghost_z).or.(periodic_bc(3))) then
    ! Get the rank of the process of coordinates [X,Y,Z+1]
    coordinates      = myid_cart
    coordinates(3)   = myid_cart(3) + 1
    call MPI_Cart_Rank(cart_comm,coordinates,target_id,ierr_MPI_rad)
    call MPI_Put(RIB%North(1,1,1), 1, RIB_XY, target_id,target_disp,1,RIB_XY,win_south,ierr_MPI_rad)
    call MPI_Get(rivb_n,           1, RIB_XY, target_id,target_disp,1,RIB_XY,win_north ,ierr_MPI_rad)
  endif

! Call Synchronisation and free memory
  call MPI_WIN_FENCE(0,win_right ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_left  ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_south ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_north ,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_bottom,ierr_MPI_rad)
  call MPI_WIN_FENCE(0,win_top   ,ierr_MPI_rad)

  call MPI_WIN_FREE(win_right ,ierr_MPI_rad)
  call MPI_WIN_FREE(win_left  ,ierr_MPI_rad)
  call MPI_WIN_FREE(win_south ,ierr_MPI_rad)
  call MPI_WIN_FREE(win_north ,ierr_MPI_rad)
  call MPI_WIN_FREE(win_bottom,ierr_MPI_rad)
  call MPI_WIN_FREE(win_top   ,ierr_MPI_rad)

! !------------------------------------------------------------------------------------------
! ! SINGLE PROCESSOR periodic BC for x direction
!   if (periodic_x==1) then
! 
! ! Left --> Right
!     rivb_r = rib%Left
! 
! ! Right --> Left
!     rivb_l = rib%Right
!
!  endif

! !------------------------------------------------------------------------------------------
! ! SINGLE PROCESSOR periodic BC for y direction
!   if (periodic_y==1) then
! 
! ! Bottom --> Top
!     rivb_t = rib%Bottom
! ! Top --> Bottom
!     rivb_b = rib%Top
! 
!   endif
! 
! !------------------------------------------------------------------------------------------
! ! SINGLE PROCESSOR periodic BC for z direction
! !
!   if (periodic_z==1) then
! 
! ! North --> South
!       rivb_s = rib%North
! ! South --> North
!       rivb_n = rib%South
! 
!    endif
! 
! 
!   PRINT*, ' EXCHANGE RIB'
!   PRINT*, 'Max rivb_r: ', maxval(rivb_r)
!   PRINT*, 'Max rivb_l: ', maxval(rivb_l)
!   PRINT*, 'Max rivb_s: ', maxval(rivb_s)
!   PRINT*, 'Max rivb_n: ', maxval(rivb_n)
!   PRINT*, 'Max rivb_t: ', maxval(rivb_t)
!   PRINT*, 'Max rivb_b: ', maxval(rivb_b)
! 
!   PRINT*,'PRINT RIB'

! Free MPI_TYPE
  call MPI_Type_free(RIB_XY,ierr_MPI_rad)
  call MPI_Type_free(RIB_YZ,ierr_MPI_rad)
  call MPI_Type_free(RIB_XZ,ierr_MPI_rad)

  return
! !----------------------------------------------------------------------------------------
  end subroutine exchange_RIVB

!========================================================================================
  subroutine create_CV
!========================================================================================
! Populate the array of controle volumes, cell_tot. cell_tot is initialize
! in subroutine allocate_radiation_arrays. Dimensions: (mx-1, my-1, mz-1)
! cell_tot is type cuboid. Recall:
!   type cuboid
!     real, dimension(2) :: x                    ! x coordinates of the cuboid
!     real, dimension(2) :: y                    ! y coordinates of the cuboid
!     real, dimension(2) :: z                    ! z coordinates of the cuboid
!     real               :: Surface              ! Surface of the control volume
!     real               :: Volume
!     real               :: Temp                 ! Temperature of the control volume (in K)
!     real               :: c_abs                ! Absorption coefficient for the control volume
!     real               :: steradiancy          ! BlackBody radiative intensity per unit of steradiant
!     real               :: Source_cv            ! Radiative source for the control volume
!     integer            :: i_x                  ! Indice of the cell in x direction
!     integer            :: i_y                  ! Indice of the cell in y direction
!     integer            :: i_z                  ! Indice of the cell in z direction
!   end type cuboid
!----------------------------------------------------------------------------------------
  
  implicit none

  real :: L_x, L_y, L_z
  
  integer :: jx, jy, jz

  cell_tot(:,:,:)%Surface     = 0.0 ! Geometrical properties
  cell_tot(:,:,:)%Volume      = 0.0 ! Geometrical properties

  cell_tot(:,:,:)%Temp        = 0.0
  cell_tot(:,:,:)%c_abs       = 0.0
  cell_tot(:,:,:)%steradiancy = 0.0
  cell_tot(:,:,:)%Source_cv   = 0.0

  do jx = 1, mx-1
    cell_tot(jx,:,:)%x(1) = xx(jx)
    cell_tot(jx,:,:)%x(2) = xx(jx+1)
  enddo

  do jy = 1, my-1
    cell_tot(:,jy,:)%y(1) = yy(jy)
    cell_tot(:,jy,:)%y(2) = yy(jy+1)
  enddo

  do jz = 1, mz-1
    cell_tot(:,:,jz)%z(1) = zz(jz)
    cell_tot(:,:,jz)%z(2) = zz(jz+1)
  enddo

  do jz = 1, mz-1
    do jy = 1, my-1
      do jx = 1, mx-1      
        cell_tot(jx,jy,jz)%Surface = surface(cell_tot(jx,jy,jz))    
        cell_tot(jx,jy,jz)%Volume  = Volume(cell_tot(jx,jy,jz))
        cell_tot(jx,jy,jz)%i_x     = jx
        cell_tot(jx,jy,jz)%i_y     = jy
        cell_tot(jx,jy,jz)%i_z     = jz
      enddo
    enddo
  enddo

  return
!----------------------------------------------------------------------------------------
  end subroutine create_CV

!========================================================================================
  subroutine Create_facets
!========================================================================================
! This function creates the boundaries surfaces of the computational domain, Boundaries
!
! This function also creates the variable Global_Cell needed for some function
! 
! There are 6 Boundaries surfaces: South, North, Right, Left, Bottom, Top
! Recall Boundaries is type face. Type face is defined as
!
! Important: This subroutine assumes that Boundary(:))%element are not allocated!
! Program will deallocate them if it is not the case
! 
!   type facet
!    type(coordinates), dimension(4) :: vertex
!    type(coordinates)               :: center
!    type(cuboid), pointer           :: cell
!    real                            :: Surface
!    integer                         :: i_dim1   ! Indices of the first dimension
!    integer                         :: i_dim2   ! Indices of the second dimension
!    type(face), pointer             :: boundary ! pointer to the boundary face if facet belongs to the a boundary
!  end type facet     
!
!   type face
!     type(facet), allocatable, dimension(:,:) :: element
!     integer                                  :: fixed_dim
!     real                                     :: Fixed_coordinate
!     real                                     :: emissivity
!     integer                                  :: plane  (1, 2, 3, 4, 5, 6)
!     integer                                  :: coordinate (1 or mz, 1 or mx, 1 or my)
!     type(coordinates), dimension(4)          :: vertex ! location of the verteces of the face
!   end type faces
!
! Populates the coordinates of the four verteces that constitute a given facet
! Compute the barycenter of the surface
! It is given by the formula:
! x_center = 0.25*(x_1+x_2+x_3+x_4)
! y_center = 0.25*(y_1+y_2+y_3+y_4)
! z_center = 0.25*(z_1+z_2+z_3+z_4)
! 
! Recall:    1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top in defined orientation
! ---------------------------------
! 1 - South : (x,y), eq: z = zz(1)
! 2 - North : (x,y), eq: z = zz(mz)
! 3 - Right : (y,z), eq: x = xx(1)
! 4 - Left  : (y,z), eq: x = xx(mx)
! 5 - Bottom: (x,z), eq: y = yy(1) 
! 6 - Top   : (x,z), eq: y = yy(my)
! 
!----------------------------------------------------------------------------------------
! 
! (i): Vertices      ______________________
!                (6)/|                    /|(7)
!                  / |  neighbor(6)      / |
!                 /  |  Plan(North)     /  |
!                /   |                 /   |
!           (5) /_____________________/    |
!              |     |                |(8) |
!              |     |          neighb|r(4)|
!              |     |            plan|top)|
! neighbor(1)  |     |                |    | neighbor(2) 
! Plan(right)  |  (2)|________________|____| Plan(left)
!              |    /                 |    /(3)
!              |   /                  |   /
!           neighbor(3)               |  /
!           Plan(bottom)              | /
!           (1)|/_____________________|/(4)
!
!                   neighbor(5)
!                   Plan(south)
!
!         (z)^  (y)
!            | /
!            |/--> (x)
! 
!----------------------------------------------------------------------------------------
 
  implicit none

 ! Coordinates of the Vertices of the computational domain
  type(coordinates), dimension(8) :: Domain_Vertices

  integer :: i_x, i_y, i_z

! Create Domain_cell: Cell that englobes the domain  

  Global_Cell%Surface     = 0.0
  Global_Cell%Volume      = 0.0
  Global_Cell%Temp        = 0.0
  Global_Cell%c_abs       = 0.0
  Global_Cell%steradiancy = 0.0
  Global_Cell%Source_cv   = 0.0
  Global_Cell%i_x = 0
  Global_Cell%i_y = 0
  Global_Cell%i_z = 0

  Global_Cell%x(1) = xx(1)
  Global_Cell%x(2) = xx(mx)

  Global_Cell%y(1) = yy(1)
  Global_Cell%y(2) = yy(my)

  Global_Cell%z(1) = zz(1)
  Global_Cell%z(2) = zz(mz)

! First allocate Boundaries()%element. Make sure they are deallocated first

  if (allocated(Boundaries(1)%element)) deallocate(Boundaries(1)%element)
  if (allocated(Boundaries(2)%element)) deallocate(Boundaries(2)%element)
  if (allocated(Boundaries(3)%element)) deallocate(Boundaries(3)%element)
  if (allocated(Boundaries(4)%element)) deallocate(Boundaries(4)%element)
  if (allocated(Boundaries(5)%element)) deallocate(Boundaries(5)%element)
  if (allocated(Boundaries(6)%element)) deallocate(Boundaries(6)%element)
  
  allocate(Boundaries(1)%element(mx-1,my-1)) ! South
  allocate(Boundaries(2)%element(mx-1,my-1)) ! North

  allocate(Boundaries(3)%element(my-1,mz-1)) ! Right
  allocate(Boundaries(4)%element(my-1,mz-1)) ! Left

  allocate(Boundaries(5)%element(mx-1,mz-1)) ! Bottom
  allocate(Boundaries(6)%element(mx-1,mz-1)) ! Top

! Create Domain_Vertices
  Domain_Vertices(1)%x = xx(1)
  Domain_Vertices(1)%y = yy(1)
  Domain_Vertices(1)%z = zz(1)

  Domain_Vertices(2)%x = xx(1)
  Domain_Vertices(2)%y = yy(my)
  Domain_Vertices(2)%z = zz(1)

  Domain_Vertices(3)%x = xx(mx)
  Domain_Vertices(3)%y = yy(my)
  Domain_Vertices(3)%z = zz(1)

  Domain_Vertices(4)%x = xx(mx)
  Domain_Vertices(4)%y = yy(1)
  Domain_Vertices(4)%z = zz(1)

  Domain_Vertices(5)%x = xx(1)
  Domain_Vertices(5)%y = yy(1)
  Domain_Vertices(5)%z = zz(mz)

  Domain_Vertices(6)%x = xx(1)
  Domain_Vertices(6)%y = yy(my)
  Domain_Vertices(6)%z = zz(mz)

  Domain_Vertices(7)%x = xx(mx)
  Domain_Vertices(7)%y = yy(my)
  Domain_Vertices(7)%z = zz(mz)

  Domain_Vertices(8)%x = xx(mx)
  Domain_Vertices(8)%y = yy(1)
  Domain_Vertices(8)%z = zz(mz)

! Assignment of fixed dimension(1 for x, 2 for y, 3 for z), 
!               Fixed_coordinate (equation of plan), emissivity
!               Vertex

! South plane

  Boundaries(1)%fixed_dim        = 3 
  Boundaries(1)%Fixed_coordinate = zz(1) 
  Boundaries(1)%emissivity       = emiss_s

  Boundaries(1)%coordinate       = 1 
  Boundaries(1)%plane            = 1 

  Boundaries(1)%vertex(1) = Domain_Vertices(1)
  Boundaries(1)%vertex(2) = Domain_Vertices(2)
  Boundaries(1)%vertex(3) = Domain_Vertices(3)
  Boundaries(1)%vertex(4) = Domain_Vertices(4)

! North plane

  Boundaries(2)%fixed_dim        = 3 
  Boundaries(2)%Fixed_coordinate = zz(mz)
  Boundaries(2)%emissivity       = emiss_n

  Boundaries(2)%coordinate       = mz 
  Boundaries(2)%plane            = 2

  Boundaries(2)%vertex(1) = Domain_Vertices(5)
  Boundaries(2)%vertex(2) = Domain_Vertices(6)
  Boundaries(2)%vertex(3) = Domain_Vertices(7)
  Boundaries(2)%vertex(4) = Domain_Vertices(8)

! Right plane

  Boundaries(3)%fixed_dim        = 1 
  Boundaries(3)%Fixed_coordinate = xx(1)
  Boundaries(3)%emissivity       = emiss_r

  Boundaries(3)%coordinate       = 1 
  Boundaries(3)%plane            = 3

  Boundaries(3)%vertex(1) = Domain_Vertices(1)
  Boundaries(3)%vertex(2) = Domain_Vertices(2)
  Boundaries(3)%vertex(3) = Domain_Vertices(6)
  Boundaries(3)%vertex(4) = Domain_Vertices(5)

! Left plane

  Boundaries(4)%fixed_dim        = 1 
  Boundaries(4)%Fixed_coordinate = xx(mx)
  Boundaries(4)%emissivity       = emiss_l

  Boundaries(4)%coordinate       = mx 
  Boundaries(4)%plane            = 4

  Boundaries(4)%vertex(1) = Domain_Vertices(4)
  Boundaries(4)%vertex(2) = Domain_Vertices(3)
  Boundaries(4)%vertex(3) = Domain_Vertices(7)
  Boundaries(4)%vertex(4) = Domain_Vertices(8)

! Bottom plane

  Boundaries(5)%fixed_dim        = 2 
  Boundaries(5)%Fixed_coordinate = yy(1)
  Boundaries(5)%emissivity       = emiss_b

  Boundaries(5)%coordinate       = 1 
  Boundaries(5)%plane            = 5

  Boundaries(5)%vertex(1) = Domain_Vertices(1)
  Boundaries(5)%vertex(2) = Domain_Vertices(4)
  Boundaries(5)%vertex(3) = Domain_Vertices(8)
  Boundaries(5)%vertex(4) = Domain_Vertices(5)

! Top plane

  Boundaries(6)%fixed_dim        = 2 
  Boundaries(6)%Fixed_coordinate = yy(my)
  Boundaries(6)%emissivity       = emiss_t

  Boundaries(6)%coordinate       = my 
  Boundaries(6)%plane            = 6

  Boundaries(6)%vertex(1) = Domain_Vertices(2)
  Boundaries(6)%vertex(2) = Domain_Vertices(3)
  Boundaries(6)%vertex(3) = Domain_Vertices(7)
  Boundaries(6)%vertex(4) = Domain_Vertices(6)

!----------------------------------------------------------------------------------------
! Assign Boundaries(:)%element(:,:)
! South and North faces

  Domain_Vertices(1:4)%z = Boundaries(1)%Fixed_coordinate
  Domain_Vertices(5:8)%z = Boundaries(2)%Fixed_coordinate

  do i_y = 1, my-1

    Domain_Vertices(1:2)%y = yy(i_y  )
    Domain_Vertices(3:4)%y = yy(i_y+1)

    Domain_Vertices(5:6)%y = yy(i_y  )
    Domain_Vertices(7:8)%y = yy(i_y+1)

    do i_x = 1, mx-1
       
      Domain_Vertices(1)%x   = xx(i_x  )
      Domain_Vertices(2:3)%x = xx(i_x+1)
      Domain_Vertices(4:5)%x = xx(i_x  )
      Domain_Vertices(6:7)%x = xx(i_x+1)
      Domain_Vertices(8)%x   = xx(i_x  )

      ! Assign vertices
      Boundaries(1)%element(i_x,i_y)%vertex(1:4) = Domain_Vertices(1:4)
      Boundaries(2)%element(i_x,i_y)%vertex(1:4) = Domain_Vertices(5:8)

      ! Compute surface of the facet 
      Boundaries(1)%element(i_x,i_y)%Surface = surface_2D(Boundaries(1)%element(i_x,i_y))
      Boundaries(2)%element(i_x,i_y)%Surface = surface_2D(Boundaries(2)%element(i_x,i_y))

      ! Compute center of the facet 
      Boundaries(1)%element(i_x,i_y)%center = compute_Center(Boundaries(1)%element(i_x,i_y))
      Boundaries(2)%element(i_x,i_y)%center = compute_Center(Boundaries(2)%element(i_x,i_y))

      ! Assign the adjacent control volume
      Boundaries(1)%element(i_x,i_y)%cell => cell_tot(i_x,i_y,1)
      Boundaries(2)%element(i_x,i_y)%cell => cell_tot(i_x,i_y,mz-1)

      ! Assign values of indices for the facets
      Boundaries(1)%element(i_x,i_y)%i_dim1 = i_x
      Boundaries(2)%element(i_x,i_y)%i_dim1 = i_x
      Boundaries(1)%element(i_x,i_y)%i_dim2 = i_y
      Boundaries(2)%element(i_x,i_y)%i_dim2 = i_y

      ! Assign values of pointer boundary
      Boundaries(1)%element(i_x,i_y)%boundary => Boundaries(1)
      Boundaries(2)%element(i_x,i_y)%boundary => Boundaries(2)

    enddo
  enddo

! Right and Left faces treatment
  Domain_Vertices(1:4)%x = Boundaries(3)%Fixed_coordinate
  Domain_Vertices(5:8)%x = Boundaries(4)%Fixed_coordinate

  do i_z = 1, mz-1

    Domain_Vertices(1:2)%z = zz(i_z)
    Domain_Vertices(3:4)%z = zz(i_z+1)

    Domain_Vertices(5:6)%z = zz(i_z)
    Domain_Vertices(7:8)%z = zz(i_z+1)

    do i_y = 1, my-1
     
      Domain_Vertices(1)%y   = yy(i_y  )
      Domain_Vertices(2:3)%y = yy(i_y+1)
      Domain_Vertices(4:5)%y = yy(i_y  )
      Domain_Vertices(6:7)%y = yy(i_y+1)
      Domain_Vertices(8)%y   = yy(i_y  )

      ! Assign vertices
      Boundaries(3)%element(i_y,i_z)%vertex(1:4) = Domain_Vertices(1:4) ! Right face
      Boundaries(4)%element(i_y,i_z)%vertex(1:4) = Domain_Vertices(5:8) ! Left  face

      ! Compute surface of the facet 
      Boundaries(3)%element(i_y,i_z)%Surface = surface_2D(Boundaries(3)%element(i_y,i_z))
      Boundaries(4)%element(i_y,i_z)%Surface = surface_2D(Boundaries(4)%element(i_y,i_z))

      ! Compute center of the facet 
      Boundaries(3)%element(i_y,i_z)%center = compute_Center(Boundaries(3)%element(i_y,i_z))
      Boundaries(4)%element(i_y,i_z)%center = compute_Center(Boundaries(4)%element(i_y,i_z))

      ! Assign the adjacent control volume
      Boundaries(3)%element(i_y,i_z)%cell => cell_tot(   1,i_y,i_z)
      Boundaries(4)%element(i_y,i_z)%cell => cell_tot(mx-1,i_y,i_z)

      ! Assign values of indices for the facets
      Boundaries(3)%element(i_y,i_z)%i_dim1 = i_y
      Boundaries(3)%element(i_y,i_z)%i_dim2 = i_z
      Boundaries(4)%element(i_y,i_z)%i_dim1 = i_y
      Boundaries(4)%element(i_y,i_z)%i_dim2 = i_z



      ! Assign values of pointer boundary
      Boundaries(3)%element(i_y,i_z)%boundary => Boundaries(3)
      Boundaries(4)%element(i_y,i_z)%boundary => Boundaries(4)

    enddo

  enddo
! Bottom and top faces
  Domain_Vertices(1:4)%y = Boundaries(5)%Fixed_coordinate
  Domain_Vertices(5:8)%y = Boundaries(6)%Fixed_coordinate

  do i_z = 1, mz-1

    Domain_Vertices(1:2)%z = zz(i_z  )
    Domain_Vertices(3:4)%z = zz(i_z+1)

    Domain_Vertices(5:6)%z = zz(i_z  )
    Domain_Vertices(7:8)%z = zz(i_z+1)

    do i_x = 1, mx-1
     
      Domain_Vertices(1)%x   = xx(i_x  )
      Domain_Vertices(2:3)%x = xx(i_x+1)
      Domain_Vertices(4:5)%x = xx(i_x  )
      Domain_Vertices(6:7)%x = xx(i_x+1)
      Domain_Vertices(8)%x   = xx(i_x  )

      ! Assign vertices
      Boundaries(5)%element(i_x,i_z)%vertex(1:4) = Domain_Vertices(1:4) ! Bottom face
      Boundaries(6)%element(i_x,i_z)%vertex(1:4) = Domain_Vertices(5:8) ! Top  face

      ! Compute surface of the facet 
      Boundaries(5)%element(i_x,i_z)%Surface = surface_2D(Boundaries(5)%element(i_x,i_z))
      Boundaries(6)%element(i_x,i_z)%Surface = surface_2D(Boundaries(6)%element(i_x,i_z))

      ! Compute center of the facet 
      Boundaries(5)%element(i_x,i_z)%center = compute_Center(Boundaries(5)%element(i_x,i_z))
      Boundaries(6)%element(i_x,i_z)%center = compute_Center(Boundaries(6)%element(i_x,i_z))

      ! Assign the adjacent control volume
      Boundaries(5)%element(i_x,i_z)%cell => cell_tot(i_x,   1,i_z)
      Boundaries(6)%element(i_x,i_z)%cell => cell_tot(i_x,my-1,i_z)

      ! Assign values of indices for the facets
      Boundaries(5)%element(i_x,i_z)%i_dim1 = i_x
      Boundaries(5)%element(i_x,i_z)%i_dim2 = i_z
      Boundaries(6)%element(i_x,i_z)%i_dim1 = i_x
      Boundaries(6)%element(i_x,i_z)%i_dim2 = i_z


      ! Assign values of pointer boundary
      Boundaries(5)%element(i_x,i_z)%boundary => Boundaries(5)
      Boundaries(6)%element(i_x,i_z)%boundary => Boundaries(6)

    enddo

  enddo
!----------------------------------------------------------------------------------------
  return
  end subroutine Create_facets


!========================================================================================
  subroutine Get_Div_qr(div_qr)
!========================================================================================
! Compute the divergence of radiative source term using the 
! cell_tot source controle volume.
!
! Variable passed in:
! None. Use global variable cell_tot$Source_cv
! 
! Variable passed out:
! div_qr, dimensions (:,:,:) Assum shape array. Units: W/m^3
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in and outer

  real, dimension(:,:,:), intent(inout) :: div_qr

! Local variables

  integer :: i, j, k
  integer :: nx, ny, nz ! Dimensions of div_qr

! Get dimensions of div_qr

  nx = ubound(div_qr,1)
  ny = ubound(div_qr,2)
  nz = ubound(div_qr,3)

!-------------------------------------------------------------------------------
! transfer radiative source needed from neighbor processors (source_cv)
! This part is more difficult than transfering temperature
! because the decomposition for the control volume is not homogeneous

!---------------------------------------------------------------------------------
! average the source of adjacent control volumes, 
! get the source at grid point (div_qr)

!----------------------
! Inner points
!----------------------

  do k = 2, nz-1
    do j = 2, ny-1
      do i = 2, nx-1
        div_qr(i,j,k)=( cell_tot(i-1,j-1,k-1)%Source_cv + cell_tot(i  ,j-1,k-1)%Source_cv +  &
                        cell_tot(i-1,j  ,k-1)%Source_cv + cell_tot(i  ,j  ,k-1)%Source_cv +  &
                        cell_tot(i-1,j-1,k  )%Source_cv + cell_tot(i  ,j-1,k  )%Source_cv +  &
                        cell_tot(i-1,j  ,k  )%Source_cv + cell_tot(i  ,j  ,k  )%Source_cv )  &
                     /( cell_tot(i-1,j-1,k-1)%Volume    + cell_tot(i  ,j-1,k-1)%Volume    +  &
                        cell_tot(i-1,j  ,k-1)%Volume    + cell_tot(i  ,j  ,k-1)%Volume    +  &
                        cell_tot(i-1,j-1,k  )%Volume    + cell_tot(i  ,j-1,k  )%Volume    +  &
                        cell_tot(i-1,j  ,k  )%Volume    + cell_tot(i  ,j  ,k  )%Volume    )          
      enddo
    enddo
  enddo


!----------------------------------------------------------------------------------------

!----------------------
! South and north plane
!----------------------

  do j = 2, ny-1
    do i = 2, nx-1

! South treatment      
!-----------------
      div_qr(i,j,1)=( cell_tot(i-1,j-1,1  )%Source_cv + cell_tot(i  ,j-1,1  )%Source_cv +     &
                      cell_tot(i-1,j  ,1  )%Source_cv + cell_tot(i  ,j  ,1  )%Source_cv )     &
                   /( cell_tot(i-1,j-1,1  )%Volume    + cell_tot(i  ,j-1,1  )%Volume    +     &
                      cell_tot(i-1,j  ,1  )%Volume    + cell_tot(i  ,j  ,1  )%Volume    )
! North Treatment
!----------------
! Recall dimension(source_cv)=(mx-1,my-1,mz-1)

      div_qr(i,j,nz)=( cell_tot(i-1,j-1,mz-1)%Source_cv + cell_tot(i  ,j-1,mz-1)%Source_cv +  &
                       cell_tot(i-1,j  ,mz-1)%Source_cv + cell_tot(i  ,j  ,mz-1)%Source_cv )  &
                    /( cell_tot(i-1,j-1,mz-1)%Volume    + cell_tot(i  ,j-1,mz-1)%Volume    +  &
                       cell_tot(i-1,j  ,mz-1)%Volume    + cell_tot(i  ,j  ,mz-1)%Volume    )
    enddo
  enddo

!----------------------
! Right and Left plane
!----------------------

  do k = 2, nz-1
    do j = 2, ny-1

! Right treatment  
!----------------
      i = 1
      div_qr(i,j,k)=( cell_tot(i  ,j-1,k-1)%Source_cv + cell_tot(i  ,j-1,k  )%Source_cv +     &
                      cell_tot(i  ,j  ,k-1)%Source_cv + cell_tot(i  ,j  ,k  )%Source_cv )     &
                   /( cell_tot(i  ,j-1,k-1)%Volume    + cell_tot(i  ,j-1,k  )%Volume    +     &
                      cell_tot(i  ,j  ,k-1)%Volume    + cell_tot(i  ,j  ,k  )%Volume    )
! Left Treatment
!----------------
      i = nx

      div_qr(i,j,k)=( cell_tot(mx-1,j-1,k-1)%Source_cv + cell_tot(mx-1,j-1,k  )%Source_cv +   &
                      cell_tot(mx-1,j  ,k-1)%Source_cv + cell_tot(mx-1,j  ,k  )%Source_cv )   &
                   /( cell_tot(mx-1,j-1,k-1)%Volume    + cell_tot(mx-1,j-1,k  )%Volume    +   &
                      cell_tot(mx-1,j  ,k-1)%Volume    + cell_tot(mx-1,j  ,k  )%Volume    )   

    enddo
  enddo

!----------------------
! Bottom and Top plane
!----------------------

  do k = 2, nz-1
    do i = 2, nx-1

! Bottom treatment  
!-----------------
      j = 1
      
      div_qr(i,j,k)=( cell_tot(i-1,j  ,k-1)%Source_cv + cell_tot(i-1,j  ,k  )%Source_cv +     &
                      cell_tot(i  ,j  ,k-1)%Source_cv + cell_tot(i  ,j  ,k  )%Source_cv )     &
                   /( cell_tot(i-1,j  ,k-1)%Volume    + cell_tot(i-1,j  ,k  )%Volume    +     &
                      cell_tot(i  ,j  ,k-1)%Volume    + cell_tot(i  ,j  ,k  )%Volume    )
! Left Treatment
!---------------
      j = ny

      div_qr(i,j,k)=( cell_tot(i-1,my-1,k-1)%Source_cv + cell_tot(i-1,my-1,k )%Source_cv +    &
                      cell_tot(i  ,my-1,k-1)%Source_cv + cell_tot(i  ,my-1,k )%Source_cv )    &
                   /( cell_tot(i-1,my-1,k-1)%Volume    + cell_tot(i-1,my-1,k )%Volume    +    &
                      cell_tot(i  ,my-1,k-1)%Volume    + cell_tot(i  ,my-1,k )%Volume    ) 

    enddo
  enddo
  
!----------------------------------------------------------------------------------------
! Edges 12 edges to do
!
! Edges along y direction
!--------------------------
  
  do j = 2,ny-1

! Edge (1,y,1)
!-------------
    div_qr(1 ,j,1)=( cell_tot(1    ,j-1, 1)%Source_cv + cell_tot(1  , j, 1)%Source_cv ) &
                  /( cell_tot(1    ,j-1, 1)%Volume    + cell_tot(1  , j, 1)%Volume    )

! Edge (nx,y,1)
!--------------
    div_qr(nx,j,1)=( cell_tot(mx-1 ,j-1, 1)%Source_cv + cell_tot(mx-1, j, 1)%Source_cv) &
                  /( cell_tot(mx-1 ,j-1, 1)%Volume    + cell_tot(mx-1, j, 1)%Volume   )

! Edge (1,y,nz)
!--------------
    
    div_qr(1,j,nz)=( cell_tot(1  ,j-1,mz-1)%Source_cv + cell_tot(1  , j,mz-1)%Source_cv) &
                  /( cell_tot(1  ,j-1,mz-1)%Volume    + cell_tot(1  , j,mz-1)%Volume   )

! Edge (nx,y,nz)
!---------------
    div_qr(nx,j,nz)=(cell_tot(mx-1,j-1,mz-1)%Source_cv + cell_tot(mx-1,j,mz-1)%Source_cv ) &
                   /(cell_tot(mx-1,j-1,mz-1)%Volume    + cell_tot(mx-1,j,mz-1)%Volume    )
  enddo

!--------------------------
! Edges along z direction
!--------------------------
  
  do k = 2,nz-1

! Edge (1,1,z)
!-------------

    div_qr( 1, 1,k)=( cell_tot(   1,  1,k-1)%Source_cv + cell_tot(  1,  1,k)%Source_cv  ) &
                   /( cell_tot(   1,  1,k-1)%Volume    + cell_tot(  1,  1,k)%Volume     )

! Edge (nx,1,z)
!--------------

    div_qr(nx, 1,k)=( cell_tot(mx-1,  1,k-1)%Source_cv + cell_tot(mx-1, 1,k)%Source_cv  ) &
                   /( cell_tot(mx-1,  1,k-1)%Volume    + cell_tot(mx-1, 1,k)%Volume     )

! Edge (1,ny,z)
!--------------
    div_qr( 1,ny,k)=( cell_tot(1  ,my-1,k-1)%Source_cv + cell_tot(1 ,my-1,k)%Source_cv  ) &
                   /( cell_tot(1  ,my-1,k-1)%Volume    + cell_tot(1 ,my-1,k)%Volume     )

! Edge (nx,ny,z)
!--------------
    div_qr(nx,ny,k)=( cell_tot(mx-1,my-1,k-1)%Source_cv + cell_tot(mx-1,my-1,k)%Source_cv ) &
                   /( cell_tot(mx-1,my-1,k-1)%Volume    + cell_tot(mx-1,my-1,k)%Volume    )

  enddo

!--------------------------
! Edges along x direction
!--------------------------
  
  do i = 2,nx-1

! Edge (x,1,1)
!-------------

    div_qr(i, 1, 1)=( cell_tot(i-1,  1,  1)%Source_cv + cell_tot(i  ,   1,  1)%Source_cv ) &
                   /( cell_tot(i-1,  1,  1)%Volume    + cell_tot(i  ,   1,  1)%Volume    )
! Edge (x,ny,1)
!--------------

    div_qr(i,ny, 1)=( cell_tot(i-1,my-1, 1)%Source_cv + cell_tot(i  ,my-1,  1)%Source_cv ) &
                   /( cell_tot(i-1,my-1, 1)%Volume    + cell_tot(i  ,my-1,  1)%Volume    )

! Edge (x,1,nz)
!--------------

    div_qr(i ,1,nz)=( cell_tot(i-1, 1,mz-1)%Source_cv + cell_tot(  i,  1,mz-1)%Source_cv ) &
                   /( cell_tot(i-1, 1,mz-1)%Volume    + cell_tot(  i,  1,mz-1)%Volume    )

! Edge (x,ny,nz)
!---------------

    div_qr(i,ny,nz)=( cell_tot(i-1,my-1,mz-1)%Source_cv + cell_tot(i,my-1,mz-1)%Source_cv ) &
                   /( cell_tot(i-1,my-1,mz-1)%Volume    + cell_tot(i,my-1,mz-1)%Volume    )

  enddo

!----------------------------------------------------------------------------------------
! Verteces

! (1,1,1)
 
  i = 1
  j = 1
  k = 1

  div_qr(i,j,k) = cell_tot(i,j,k)%Source_cv/cell_tot(i,j,k)%Volume
   
! (1,ny,1)
  
  i = 1
  j = ny
  k = 1

  div_qr(i,j,k) = cell_tot(i,my-1,k)%Source_cv/cell_tot(i,my-1,k)%Volume

! (nx,1,1)
  
  i = nx
  j = 1
  k = 1

  div_qr(i,j,k) = cell_tot(mx-1,j,k)%Source_cv/cell_tot(mx-1,j,k)%Volume

! (nx,ny,1)

  i = nx
  j = ny
  k = 1

  div_qr(i,j,k) = cell_tot(mx-1,my-1,k)%Source_cv/cell_tot(mx-1,my-1,k)%Volume

! (1,1,nz)
  
  i = 1
  j = 1
  k = nz

  div_qr(i,j,k) = cell_tot(i,j,mz-1)%Source_cv/cell_tot(i,j,mz-1)%Volume
   
! (1,ny,nz)

  j = ny

  div_qr(i,j,k) = cell_tot(1,my-1,mz-1)%Source_cv/cell_tot(1,my-1,mz-1)%Volume

! (nx,1,nz)
  
  i = nx
  j = 1
  k = nz

  div_qr(i,j,k) = cell_tot(mx-1,j,mz-1)%Source_cv/cell_tot(mx-1,j,mz-1)%Volume

! (nx,ny,nz)

  i = nx
  j = ny
  k = nz

  div_qr(i,j,k) = cell_tot(mx-1,my-1,mz-1)%Source_cv/cell_tot(mx-1,my-1,mz-1)%Volume

! END OF TRANSFER RAD SOURCE, TRANSFER VOLUME-POINT
!--------end of distributing source to grid point -------------------------------------
  return
  end subroutine Get_Div_qr

!========================================================================================
  function Get_Points_Index(i,j,k,origin) result(j_point)
!========================================================================================

! Recall : 
! 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top 
  implicit none
  integer, intent(in) :: i, j, k, origin
  integer :: j_point

  select case(origin)

    case(1) ! South
      j_point = (i-1)*(my-1)+j
    case(2) ! North
      j_point = (i-1)*(my-1)+j+(mx-1)*(my-1)
    case(3) ! Right
      j_point = (j-1)*(mz-1)+k+2*(mx-1)*(my-1)
    case(4) ! Left
      j_point = (j-1)*(mz-1)+k+2*(mx-1)*(my-1)+(my-1)*(mz-1)
    case(5) ! Bottom
      j_point = (k-1)*(mx-1)+i+2*(mx-1)*(my-1)+2*(my-1)*(mz-1)
    case(6) ! Top
      j_point = (k-1)*(mx-1)+i+2*(mx-1)*(my-1)+2*(my-1)*(mz-1)+(mz-1)*(mx-1)
    case default
      j_point = -1
      print*, 'Origin plane out of range'
  end select 

  return
!----------------------------------------------------------------------------------------
  end function Get_Points_Index

!========================================================================================
  subroutine Find_Associated_Ray(DTM_rays)
!========================================================================================
! The purpose of this subroutine is to populate the "iray2" attribute 
! of variable Point_origin (type location)
! It corresponds to the ray that has the closest direction vector to the ray treated and 
! located at the location of origin. This is important for virtual boundaries and 
! periodic boundary conditions
! 
! Variables passed in & outer
! DTM_ray :: type(beam), dimension(:,:): First dimension:  overall number of facets, (=n_origin(1))
!                                        Second dimension: numbers of rays traced from each
!                                        facets (=n_origin(2))
!
!----------------------------------------------------------------------------------------
!  Recall
!========
!        top:    yy(my)
!        bottom: yy(1)
!
!        left:   xx(mx)
!        right:  xx(1)
!
!        north:  zz(mz)
!        south:  zz(1)
!        1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top 
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in and out
  type(beam), intent(INOUT), dimension(:,:), target :: DTM_rays

! Local variables
  type(beam), pointer :: single_ray

  real, pointer :: Phi_theta(:)

  real, dimension(3) :: vec_dir1
  real, dimension(3) :: vec_dir2

  real :: theta1, theta2, phi1, phi2

  integer, dimension(2) :: n_origin
  integer, dimension(4) :: Index_facet
  
  integer :: i_phi, i_theta
  integer :: i_phi2, i_theta2
  integer :: i, k
  integer :: i_ray, i_ray2 
  integer :: j_point 
  integer :: edge, coordin_origin

! Get dimensions of Point_origin

  n_origin = SHAPE(DTM_rays)
  
! recall: if DTM_Ray is global variable ray, then we have
! n_origin(1) = n_rays
! n_origin(2) = n_points

!------------------------------------------------------------------------------
! Loop over the elements of the variable DTM_rays

  do j_point = 1, n_origin(2)
    do i_ray = 1, n_origin(1)

      single_ray => DTM_rays(i_ray,j_point)

      edge    = single_ray%edge_final ! Boundary plane of point j_point

      i_phi   = single_ray%iphi
      i_theta = single_ray%itheta

!------------------------------------------------------------------------------
! Get the indices of the facet "origin" for ray "single_ray"

      Index_facet = Get_Facet_Coordinates(single_ray%Facet_origin)

! Determination of the reflective plane for virtual boundary conditions
! This is needed to find the associated ray (iray2) at location of virtual boundary

      select case(Index_facet(4))
        case(1)  ! South plane
          coordin_origin = 2
        case(2)  ! North plane
          coordin_origin = 1
        case(3)  ! Right plane
          coordin_origin = 4
        case(4)  ! Left plane
          coordin_origin = 3
        case(5)  ! Bottom plane
          coordin_origin = 6
        case(6)  ! Top plane
          coordin_origin = 5
      end select


!------------------------------------------------------------------------------
! Find the associated ray (iray2) at the origin for virtual boundary and periodic
! boundary conditions. This is done at initialization.

! Initialization i_theta2
      i_theta2 = 0
      i_phi2   = 0

      theta1 = theta(i_theta)
      phi1   = phi(i_theta)%azimuth(i_phi)

! Compute direction vector for ith ray. The coordinates are expressed 
! in the base associated with plane "edge"

      vec_dir1 = 0.0

      vec_dir1(1) = cos(phi1)*sin(theta1)
      vec_dir1(2) = sin(phi1)*sin(theta1)
      vec_dir1(3) = cos(theta1)

      vec_dir2 = 0.0

! Perfom change of base. Express vec_dir1 in base associated with plane "coordin_origin"
  
      do k = 1,3
        vec_dir2(k)=DOT_PRODUCT(Change_plane(edge,coordin_origin,k,1:3),vec_dir1(1:3))
      end do 

! atan2(y,x) gives angle phi For any real arguments x and y not both equal to zero,
! atan2(y, x) is the angle in radians between the positive x-axis of a plane and the
! point given by the coordinates (x, y) on it. The angle is positive for counter-clockwise
! angles (upper half-plane, y > 0), and negative for clockwise angles
!  (lower half-plane, y < 0). 
! Phi2 is between 0 and from 0 to 2 pi (excluded)

      phi2   = atan2(vec_dir2(2),vec_dir2(1)) 
      theta2 = acos(vec_dir2(3))

! Test if phi2 is negative. If yes then do 2pi+phi2
      if (phi2.lt.0.0) phi2 = 2.0*pi + phi2
      if (phi2.ge.2.0*pi) phi2 = phi2 - 2.0*pi

!------------------------------------------------------------------------------
! Find nearest theta and phi
!
! Phi2 Find index of the nearest equivalent phi
! if phi2 <= phi(1), then assume phi2 = phi(1)
! elseif phi(n_phi) <= phi2, assume phi2 = phi(n_phi)
! elseif case (phi(1) <= phi2 <= phi(n_phi)
! if (0.5*(phi(i)+phi(i+1)).lt.phi2) then phi2 is closer to phi(i_phi2+1)
!  i_phi2 = i+1
! else
!  i_phi2 = i

! Theta2 Find index of the nearest equivalent theta
FindTheta: if (theta2.le.theta(1)) then            !  theta2 <= theta(1), then assume theta2 = theta(1) 
        i_theta2 = 1
      elseif (theta2.ge.theta(n_theta)) then 	!  theta(n_theta) <= theta2, assume theta2 = theta(n_theta) 
        i_theta2 = n_theta
      else                                        ! case (theta(1) <= theta2 <= theta(n_theta)
        do i = 1, n_theta-1
          if ((theta(i)-theta2)*(theta(i+1)-theta2).le.0.0) then  ! theta(i_theta2)<=theta2<=theta(i_theta2+1)
            if (0.5*(theta(i)+theta(i+1)).lt.theta2) then ! theta2 is closer to theta(i_theta2+1)
              i_theta2 = i+1
            else
              i_theta2 = i
            endif
            exit  ! exit loop, nearest theta found
          endif   ! theta2< theta(i_theta2) or theta(i_theta2+1) < theta2
          ! BLANK KEEP LOOPING    
        enddo
      endif FindTheta

      Phi_theta => phi(i_theta2)%azimuth

FindPhi:if (phi2.le.Phi_theta(1)) then    			!  phi2 <= Phi_theta(1), then assume phi2 = Phi_theta(1) 
        i_phi2 = 1
      elseif (phi2.ge.Phi_theta(n_phi(i_theta2))) then  !  Phi_theta(n_phi(i_theta2)) <= phi2, assume phi2 = Phi_theta(n_phi(i_theta2)) 
        i_phi2 = n_phi(i_theta2)
      else                              ! case (Phi_theta(1) <= phi2 <= Phi_theta(n_phi(i_theta2))
        do i = 1, n_phi(i_theta2)-1
          if ((Phi_theta(i)-phi2)*(Phi_theta(i+1)-phi2).le.0.0) then  ! Phi_theta(i_phi2)<=phi2<=Phi_theta(i_phi2+1)
            if (0.5*(Phi_theta(i)+Phi_theta(i+1)).lt.phi2) then ! phi2 is closer to Phi_theta(i_phi2+1)
              i_phi2 = i+1
            else
              i_phi2 = i
            endif
            exit ! exit loop, nearest phi found
          endif  ! phi2< Phi_theta(i_phi2) or Phi_theta(i_phi2+1) < phi2
        ! BLANK KEEP LOOPING    
        enddo
      endif FindPhi

!------------------------------------------------------------------------------
! Compute ray indice associated with i_theta2 and i_phi2

      i_ray2 = sum(n_phi(1:i_theta2-1))+i_phi2

!------------------------------------------------------------------------------
! Save the associated ray 

      single_ray%iray2 = i_ray2

!------------------------------------------------------------------------------
! Deallocate pointers
      nullify(Phi_theta)
      nullify(single_ray)

    enddo
  enddo

  return
!----------------------------------------------------------------------------------------
  end subroutine Find_Associated_Ray

!========================================================================================
  subroutine Init_MPI_Rad
!========================================================================================
! This subroutine initialize the MPI variables needed for the execution of the Radiation
! module
  use param_m,     only : periodic_x,periodic_y,periodic_z
  use topology_m,  only : myid, gcomm, ierr, xpes, ypes, zpes

  include 'mpif.h'

! Local variables

  integer, parameter         :: ndim = 3 ! Number of dimensions
  integer, dimension(2*ndim) :: neighbor
  integer, dimension(ndim)   :: direction

  integer, dimension(2*ndim) :: Cart_2_orientation ! Vector for translation from MPI 
                                                   ! Cartesian topology to 
                                                   ! defined orientation:
                                                   ! South, North, etc...

  real :: buffer

  integer :: shift
  integer :: i
!----------------------------------------------------------------------------------------
! Define Cart_2_orientation
!
! Recall:    1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top in defined orientation
! Translated:5      , 6      , 1      , 2     , 3,        4 

  Cart_2_orientation = [3,4,5,6,1,2]
!----------------------------------------------------------------------------------------
! create Cartesian coordinates for use with Radiation module

  periodic_bc = .FALSE.

  if (periodic_x == 1) periodic_bc(1) = .TRUE.
  if (periodic_y == 1) periodic_bc(2) = .TRUE.
  if (periodic_z == 1) periodic_bc(3) = .TRUE.

  cart_dim(1) = xpes
  cart_dim(2) = ypes
  cart_dim(3) = zpes

  ierr_MPI_rad = ierr

  call MPI_Cart_Create(gcomm,ndim,cart_dim,periodic_bc,.false.,cart_comm,ierr_MPI_rad) 

! Get CARTESIAN COORDINATES in communicator cart_dim
! Store it in myid_cart

  call MPI_Cart_coords(cart_comm,myid,ndim,myid_cart,ierr_MPI_rad) 

! Define ghost_x, ghost_y, ghost_z using MPI_CART_SHIFT
! Get the rank of the neighboring process in the chosen direction  increased by ones
  shift = 1
  direction(1) = 0 ! Direction along x
  direction(2) = 1 ! Direction along y
  direction(3) = 2 ! Direction along z
  
  do i = 1, ndim
    call MPI_Cart_Shift(cart_comm, direction(i), shift, neighbor(2*(i-1)+1), neighbor(2*(i-1)+2), ierr_MPI_rad)
  enddo

! Test if ghost plane needed. Remember that virtual cartesian topology implicitly accounts 
! for periodic boundary conditions

! Initialize variable is_vb
  do i = 1, 6
      is_vb(Cart_2_orientation(i)) = (neighbor(i).ge.0)
  enddo

! Test if rank of process with Cartesian coordinates [X+1,Y,Z] exists or note
! if it does not, then it means we have an open boundary, no need for ghost planes
! Do not test the x-most processor
  if (myid_cart(1).ge.(cart_dim(1)-1)) then 
    ghost_x = .false.
  else
    ghost_x = .true.
  endif

! Perform same test but in y direction
! Do not test the y-most processor
  if (myid_cart(2).ge.(cart_dim(2)-1)) then 
    ghost_y = .false.
  else
    ghost_y = .true.
  endif

! Perform same test but in z direction
! Do not test the z-most processor
  if (myid_cart(3).ge.(cart_dim(3)-1)) then 
    ghost_z = .false.
  else
    ghost_z = .true.
  endif

! Get size of Variable element type  
!  call MPI_SIZEOF(buffer,sizereal,ierr_MPI_rad)
  select Case (KIND(buffer))
    case(4)
      call MPI_Type_Size(MPI_REAL,size_type_real,ierr_MPI_rad)
      call MPI_Type_Contiguous(1,MPI_REAL,TypeReal,ierr_MPI_rad)
      call MPI_Type_commit(TypeReal,ierr_MPI_rad)
    case(8)
      call MPI_Type_Size(MPI_DOUBLE,size_type_real,ierr_MPI_rad)
      call MPI_Type_Contiguous(1,MPI_DOUBLE,TypeReal,ierr_MPI_rad)
      call MPI_Type_commit(TypeReal,ierr_MPI_rad)
    case default
      call MPI_ABORT(gcomm,ierr_MPI_rad) ! REAL FORMAT NOT SUPPORTED BY MPI. ABORT
  end select  

  return
!----------------------------------------------------------------------------------------
  end subroutine Init_MPI_Rad

!========================================================================================
  subroutine write_run_attribute(group_id,io)
!========================================================================================
! This subroutine writes the metadata
! The data written follows the format:
!----------------------------------------------------------------------------------------  
! Create group attribute. Save metadata.
! Case title : run_title (CHARACTER)
! Date of writing : value(1) (year), value(2) (month), value(3) (day)
! Time of writing : value(5) (hour), value(6) (min), value(7) (sec) 'GMT ZONE(1:2)hZONE(3:4)min
! # proc in x     : cart_dim(1) (INTEGER) - radiation_m
! # proc in y     : cart_dim(1) (INTEGER) - radiation_m
! # proc in z     : cart_dim(1) (INTEGER) - radiation_m
! Dimension                   : numdim     (INTEGER) - param_m
! Periodic in x               : periodic_x (INTEGER) - param_m
! Periodic in y               : periodic_y (INTEGER) - param_m
! Periodic in z               : periodic_z (INTEGER) - param_m
! number of species (without N2) : nsc  (INTEGER) - param_m
!----------------------------------------------------------------------------------------
  use param_m,  only : numdim, periodic_x, periodic_y,  periodic_z, nsc
  use HDF5

  implicit none

! Variable passed IN
  integer(HID_T), intent(INOUT) :: group_id      ! group identifier
  integer, intent(IN)           :: io 
! Local variables

  real                  :: atr_real_data
  integer, dimension(8) :: value
  character(LEN=100)    :: attribute_name  ! Attribute name
  character(LEN=100)    :: atr_data        ! Attribute data

  character(LEN=10)  date, time_char, zone, login


  integer :: npx, npy, npz

  npx = cart_dim(1)
  npy = cart_dim(2)
  npz = cart_dim(3)

! Get login, date, and time for records as File attribute

!  call getlog(login)
  call date_and_time(date,time_char,zone,value)
 
  if (login=='') login='N/A'

!----------------------------------------------------------------------------------------
  attribute_name = 'USER'
  atr_data       = login
  call write_STRING_attribute(group_id,attribute_name, atr_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'DATE'
  atr_data       =  date(5:6)//'/'//date(7:8)//'/'//date(1:4)
  call write_STRING_attribute(group_id,attribute_name, atr_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'TIME'
  atr_data       =  time_char(1:2) //'H'//time_char(3:4)//'MIN'//' GMT'//zone(1:3)
  call write_STRING_attribute(group_id,attribute_name, atr_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'N-PROC IN X'
  atr_real_data  =  float(npx)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'N-PROC IN Y'
  atr_real_data  =  float(npy)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'N-PROC IN Z'
  atr_real_data  =  float(npz)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'N-DIMENSIONS'
  atr_real_data  =  float(numdim)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'X-PERIODIC'
  atr_real_data  =  float(periodic_x)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'Y-PERIODIC'
  atr_real_data  =  float(periodic_y)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'Z-PERIODIC'
  atr_real_data  =  float(periodic_z)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!----------------------------------------------------------------------------------------
  attribute_name = 'TOTAL SPECIES'
  atr_real_data  =  float(nsc+1)
  call write_REAL_attribute(group_id,attribute_name, atr_real_data,io)

  return
!----------------------------------------------------------------------------------------
  end subroutine write_run_attribute

!========================================================================================
  subroutine read_run_attribute(group_id,io)
!========================================================================================
! This subroutine reads the metadata needed for restart
! The data read follow the format:
!----------------------------------------------------------------------------------------  
! Read following attributes for restart.
! # proc in x     : npx (INTEGER) - param_m
! # proc in y     : npy (INTEGER) - param_m
! # proc in z     : npz (INTEGER) - param_m
!----------------------------------------------------------------------------------------
  use param_m,  only : numdim, periodic_x, periodic_y,  periodic_z, nsc
  use HDF5

  implicit none

! Variable passed IN
  integer(HID_T), intent(INOUT) :: group_id      ! group identifier
  integer, intent(IN)           :: io 
! Local variables

  real                  :: atr_real_data
  integer, dimension(8) :: value
  character(LEN=100)    :: attribute_name  ! Attribute name
  character(LEN=100)    :: atr_data        ! Attribute data

  integer :: npx, npy, npz

  npx = cart_dim(1)
  npy = cart_dim(2)
  npz = cart_dim(3)


!----------------------------------------------------------------------------------------
  attribute_name = 'N-PROC IN X'
  call read_REAL_attribute(group_id,attribute_name, atr_real_data,io)
  
  if (int(atr_real_data)/=npx) then
    write(io,*) 'The number of processors in X has been changed!'
    write(io,*) 'From ', int(atr_real_data), ' to ', npx 
  endif
!----------------------------------------------------------------------------------------
  attribute_name = 'N-PROC IN Y'
  call read_REAL_attribute(group_id,attribute_name, atr_real_data,io)

  if (int(atr_real_data)/=npy) then
    write(io,*) 'The number of processors in Y has been changed!'
    write(io,*) 'From ', int(atr_real_data), ' to ', npy 
  endif
!----------------------------------------------------------------------------------------
  attribute_name = 'N-PROC IN Z'
  call read_REAL_attribute(group_id,attribute_name, atr_real_data,io)

  if (int(atr_real_data)/=npz) then
    write(io,*) 'The number of processors in Z has been changed!'
    write(io,*) 'From ', int(atr_real_data), ' to ', npz 
  endif

  return
!----------------------------------------------------------------------------------------
  end subroutine read_run_attribute


!========================================================================================
  subroutine write_STRING_attribute(group_id,attribute_name, atr_data,io)
!========================================================================================
! This subroutine is a private subroutine. Not intended for user call.
! Create an attribute <attribute_name> of type STRING, with values <atr_data> 
! in group <group_id>. In case of error, write in io.
! Returns nothing

  use HDF5
  use topology_m,  only : hdf5_error

  implicit none
! Variables passed in
  integer(HID_T), intent(INOUT) :: group_id        ! group identifier
  character(LEN=100),intent(IN) :: attribute_name  ! Attribute name
  character(LEN=100),intent(IN) :: atr_data        ! Attribute data
  integer, intent(IN)           :: io 

! Local variables

  integer(HID_T)  :: attr_id       ! Attribute identifier
  integer(HID_T)  :: aspace_id     ! Attribute Dataspace identifier
  integer(HID_T)  :: atype_id      ! Attribute Dataspace identifier
  integer(HSIZE_T), dimension(1) :: adims = (/1/) ! Attribute dimension
  integer         :: arank = 1                    ! Attribure rank
  integer(SIZE_T) :: attrlen       ! Length of the attribute string


  attrlen = LEN_TRIM(atr_data)

  call h5tcopy_f(H5T_NATIVE_CHARACTER, atype_id, hdf5_error)
  call h5tset_size_f(atype_id, attrlen, hdf5_error)

  call h5screate_simple_f(arank, adims, aspace_id, hdf5_error)

  call h5acreate_f(group_id, trim(attribute_name), atype_id, aspace_id, &
                      attr_id, hdf5_error)

  call h5awrite_f(attr_id, atype_id,trim(atr_data), adims, hdf5_error)
  call h5aclose_f(attr_id, hdf5_error)

  return
!----------------------------------------------------------------------------------------
  end subroutine write_STRING_attribute
  
!========================================================================================
  subroutine write_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!========================================================================================
! This subroutine is a private subroutine. Not intended for user call.
! Create an attribute <attribute_name> of type STRING, with values <atr_data> 
! in group <group_id>. In case of error, write in io.
! Returns nothing

  use HDF5
  use topology_m,  only : hdf5_error, myid

  implicit none
! Variables passed in
  integer(HID_T), intent(INOUT) :: group_id        ! group identifier
  character(LEN=100),intent(IN) :: attribute_name  ! Attribute name
  real,intent(IN)               :: atr_real_data   ! Attribute data
  integer, intent(IN)           :: io 

! Local variables
  integer(HID_T) :: write_real_type   ! Type of data to write (REAL, DOUBLE, etc...)
  integer(HID_T)  :: attr_id       ! Attribute identifier
  integer(HID_T)  :: aspace_id     ! Attribute Dataspace identifier
  integer(HSIZE_T), dimension(1) :: adims = (/1/) ! Attribute dimension
  integer         :: arank = 1                    ! Attribure rank
  integer(SIZE_T) :: attrlen       ! Length of the attribute string

  if (KIND(atr_real_data) == 4) then      ! Get type of field, Could be real (Kind = 4) 
    write_real_type = H5T_NATIVE_REAL  ! or double (Kind = 8). Need for HDF5 subroutines
  else if(KIND(atr_real_data) == 8) then
    write_real_type = H5T_NATIVE_DOUBLE
  else
    if (myid == 0) then
      write(io,*) 'The datatype to be written with HDF5 subroutine should be either REAL or DOUBLE.'
      write(io,*) 'ERROR in write_Field_dataset. Wrong datatype. Skip writting procedure.'         
      return
    endif
  endif


  call h5screate_simple_f(arank, adims, aspace_id, hdf5_error)

  call h5acreate_f(group_id, trim(attribute_name), write_real_type, aspace_id, &
                      attr_id, hdf5_error)

  call h5awrite_f(attr_id, write_real_type,atr_real_data, adims, hdf5_error)
  call h5aclose_f(attr_id, hdf5_error)

  return
!----------------------------------------------------------------------------------------
  end subroutine write_REAL_attribute

!========================================================================================
  subroutine read_REAL_attribute(group_id,attribute_name, atr_real_data,io)
!========================================================================================
! This subroutine is a private subroutine. Not intended for user call.
! Create an attribute <attribute_name> of type STRING, with values <atr_data> 
! in group <group_id>. In case of error, write in io.
! Returns nothing

  use HDF5
  use topology_m,  only : hdf5_error, myid

  implicit none
! Variables passed in
  integer(HID_T), intent(INOUT) :: group_id        ! group identifier
  character(LEN=100),intent(IN) :: attribute_name  ! Attribute name
  real,intent(OUT)              :: atr_real_data   ! Attribute data
  integer, intent(IN)           :: io 

! Local variables
  integer(HID_T)  :: read_real_type ! Type of data to read (REAL, DOUBLE, etc...)
  integer(HID_T)  :: attr_id        ! Attribute identifier
  integer(HSIZE_T), dimension(1) :: adims = (/1/) ! Attribute dimension
  logical         :: attr_exists    ! if TRUE, then the attribute exist

  if (KIND(atr_real_data) == 4) then      ! Get type of field, Could be real (Kind = 4) 
    read_real_type = H5T_NATIVE_REAL  ! or double (Kind = 8). Need for HDF5 subroutines
  else if(KIND(atr_real_data) == 8) then
    read_real_type = H5T_NATIVE_DOUBLE
  else
    if (myid == 0) then
      write(io,*) 'The datatype to be read with HDF5 subroutine should be either REAL or DOUBLE.'
      write(io,*) 'ERROR in write_Field_dataset. Wrong datatype. Skip writting procedure.'         
      return
    endif
  endif

  call h5aexists_f(group_id, trim(attribute_name),attr_exists,hdf5_error)
  
  if (attr_exists) then 
! Attribute <attribute_name> exists. Open and get its value.

    call h5aopen_f(group_id, trim(attribute_name), attr_id, hdf5_error)
    call h5aread_f(attr_id, read_real_type, atr_real_data, adims, hdf5_error)
    call h5aclose_f(attr_id, hdf5_error)
  else
    if (myid == 0) then 
      write(io,*) 'The attribute named '//trim(attribute_name) // 'does NOT exists.'
    endif
  endif
  
  return
!----------------------------------------------------------------------------------------
  end subroutine read_REAL_attribute

!========================================================================================
  subroutine write_File_dataset(file_id,field,field_dim,field_name,unit,conversion,io)
!========================================================================================
  use topology_m, only : myid, hdf5_error
  use param_m,    only : nx_g, ny_g, nz_g

  use HDF5

  implicit none

  integer(HID_T), intent(IN)            :: file_id      ! File identifier
  integer, dimension(3), intent(IN)     :: field_dim    ! Vector of fiel dimension
  real, dimension(field_dim(1),field_dim(2),field_dim(3)), intent(IN) :: field ! Data to be stored
  character(LEN=100), intent(IN)        :: field_name   ! Name of the data_set
  character(LEN=100), intent(IN)        :: unit         ! String identifying the unit
  real, intent(IN)                      :: conversion   ! Constant to multiply field with
  integer, intent(INOUT)                :: io

!-- local 

!--HDF5 Variables
  integer(HID_T) :: dset_id       ! Dataset identifier 
  integer(HID_T) :: filespace     ! Dataspace identifier in file 
  integer(HID_T) :: memspace      ! Dataspace identifier in memory
  integer(HID_T) :: plist_id      ! Property list identifier
  integer(HID_T) :: write_type    ! Type of data to write (REAL, DOUBLE, etc...)


! HDF5 variables for attribute
  character(LEN=100),parameter :: attribute_name1 = "UNITS" ! Attribute name
  character(LEN=100),parameter :: attribute_name2 = "CONV_COEFF" ! Attribute name

  integer(HID_T)  :: attr_id       ! Attribute identifier
  integer(HID_T)  :: aspace_id     ! Attribute Dataspace identifier
  integer(HID_T)  :: atype_id      ! Attribute Dataspace identifier
  integer(HSIZE_T), dimension(1) :: adims = (/1/) ! Attribute dimension
  integer         :: arank = 1                    ! Attribure rank
  integer(SIZE_T) :: attrlen       ! Length of the attribute string

  integer(HSIZE_T), allocatable :: dimsf(:) ! Dataset dimensions
                                            ! in the file.
  integer(HSIZE_T),  allocatable :: chunk_dims(:) ! Chunks dimensions
  integer(HSIZE_T),  allocatable :: count_ele(:)  
  integer(HSSIZE_T), allocatable :: offset(:) 
  integer(HSIZE_T) , allocatable :: stride(:)
  integer(HSIZE_T) , allocatable :: block(:)

  real, allocatable, target       :: reshape_field(:,:,:)
  real, pointer, dimension(:)     :: point_1D
  real, pointer, dimension(:,:)   :: point_2D
  real, pointer, dimension(:,:,:) :: point_3D

  integer :: rank
  integer :: i ! dummy index

  integer :: xid, yid, zid

  xid = myid_cart(1)
  yid = myid_cart(2)
  zid = myid_cart(3)

  ! Define rank of the array field

  rank = 0

  do i = 1, 3
    if (field_dim(i)>1) rank = rank + 1
  enddo

  if (rank == 0) then
    rank = 1
  endif

!----------------------------------------------------------------------------------------
  if (KIND(field) == 4) then      ! Get type of field, Could be real (Kind = 4) 
    write_type = H5T_NATIVE_REAL  ! or double (Kind = 8). Need for HDF5 subroutines
  else if(KIND(field) == 8) then
    write_type = H5T_NATIVE_DOUBLE
  else
    if (myid == 0) then
      write(io,*) 'The datatype to be written with HDF5 subroutine should be either REAL or DOUBLE.'
      write(io,*) 'ERROR in write_Field_dataset. Wrong datatype. Skip writting procedure.'         
      return
    endif
  endif

  ! Allocate variables needed for I/O
  allocate(dimsf(rank))
  allocate(chunk_dims(rank))
  allocate(count_ele(rank))
  allocate(offset(rank))
  allocate(stride(rank))
  allocate(block(rank))

  allocate(reshape_field(field_dim(3),field_dim(2),field_dim(1)))
  reshape_field = reshape(field,(/field_dim(3),field_dim(2),field_dim(1)/),ORDER=(/3,2,1/))       

  ! Define shape of dataset in file and memory
  ! Set up for 2D only, Need addition of extra case for 3D, i.e. case when 
  ! field_dims(3)/=1 when all the other are set to 1

  Select case(rank)
    case(3)

      dimsf      = (/nz_g, ny_g, nx_g/)
      chunk_dims = (/field_dim(3),field_dim(2),field_dim(1)/)
      point_3D   => reshape_field

      stride(1) = 1 
      stride(2) = 1
      stride(3) = 1 

      count_ele(1)  =  1 
      count_ele(2)  =  1
      count_ele(3)  =  1

      block(1)  = chunk_dims(1)
      block(2)  = chunk_dims(2)
      block(3)  = chunk_dims(3)

      offset(1) = zid*chunk_dims(1)
      offset(2) = yid*chunk_dims(2)
      offset(3) = xid*chunk_dims(3)
    case(2)  
      dimsf      = (/ny_g, nx_g/)
      chunk_dims = (/field_dim(2),field_dim(1)/)
      point_2D   => reshape_field(1,:,:)
      stride(1) = 1 
      stride(2) = 1 
  
      count_ele(1)  =  1 
      count_ele(2)  =  1 
      
      block(1)  = chunk_dims(1)
      block(2)  = chunk_dims(2)

      offset(1) = yid*chunk_dims(1)
      offset(2) = xid*chunk_dims(2)

    case(1)
      if (field_dim(1)/=1) then
        point_1D  => reshape_field(1,1,:)
        dimsf        = (/nx_g/)
        chunk_dims   = (/field_dim(1)/)
        count_ele(1) = 1
        offset(1)    = xid*chunk_dims(1)
        stride(1)    = 1
        block(1)     = chunk_dims(1)
      elseif (field_dim(2)/=1) then
        point_1D  => reshape_field(1,:,1)
        dimsf        = (/ny_g/)
        chunk_dims   = (/field_dim(2)/)
        count_ele(1) = 1
        offset(1)    = yid*chunk_dims(1)
        stride(1)    = 1
        block(1)     = chunk_dims(1)
      else
        point_1D  => reshape_field(:,1,1)
        dimsf        = (/nz_g/)
        chunk_dims   = (/field_dim(3)/)
        count_ele(1) = 1
        offset(1)    = zid*chunk_dims(1)
        stride(1)    = 1
        block(1)     = chunk_dims(1)
      endif
      
  end select

  ! allocate arrays with corresponding number of dimensions

  call h5screate_simple_f(rank, dimsf, filespace, hdf5_error)
  call h5screate_simple_f(rank, chunk_dims, memspace, hdf5_error)
  !
  ! Create chunked dataset.
  !
  call h5pcreate_f(H5P_DATASET_CREATE_F, plist_id, hdf5_error)
  call h5pset_chunk_f(plist_id, rank, chunk_dims, hdf5_error)

  call h5dcreate_f(file_id, trim(field_name), write_type, filespace, &
                      dset_id, hdf5_error, plist_id)

  ! Close data space created to free memory
  !
  call h5sclose_f(filespace, hdf5_error)

  call h5dget_space_f(dset_id, filespace, hdf5_error)
  call h5sselect_hyperslab_f (filespace, H5S_SELECT_SET_F, offset, count_ele, hdf5_error, &
                                 stride, block)
  !
  ! Create property list for collective dataset write
  !
  call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdf5_error) 
  call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdf5_error)     
  !
  ! Write the dataset collectively. 
  !
  select case (rank)
    case(3)
      call h5dwrite_f(dset_id, write_type, point_3D*conversion, dimsf, hdf5_error, &
                    file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)
    case(2)
      call h5dwrite_f(dset_id, write_type, point_2D*conversion, dimsf, hdf5_error, &
                    file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)
    case(1)
      call h5dwrite_f(dset_id, write_type, point_1D*conversion, dimsf, hdf5_error, &
                    file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)
  end select


  ! Create the attributes for dataset. Print the unit and the dimensionalizing 
  ! coefficient

  call write_STRING_attribute(dset_id,attribute_name1,unit,io)
  call write_REAL_attribute(dset_id,attribute_name2,conversion,io)

  ! Close dataspaces.
  call h5sclose_f(filespace, hdf5_error)
  call h5sclose_f(memspace, hdf5_error)
  ! Close the dataset.
  call h5dclose_f(dset_id, hdf5_error)
  ! Close the property list.
  call h5pclose_f(plist_id, hdf5_error)

  select case (rank)
    case(3)
      nullify(point_3D)
    case(2)
      nullify(point_2D)
    case(1)
      nullify(point_1D)
   end select
  
  deallocate(dimsf)
  deallocate(chunk_dims)
  deallocate(count_ele)
  deallocate(offset)
  deallocate(stride)
  deallocate(block)

  return
!----------------------------------------------------------------------------------------
  end subroutine write_File_dataset

!========================================================================================
  subroutine read_File_dataset(file_id,field,field_dim,field_name,unit,conversion,io)
!========================================================================================
! READ FILE IN HFD5 FORMAT FILES
! INPUT:
! file_id    :
! field_dim  :
! field_name :
! conversion :
! io         :
! OUTPUT     
! field      :
  use topology_m, only : myid, hdf5_error
  use param_m,    only : nx_g, ny_g, nz_g

  use HDF5

  implicit none

  integer(HID_T), intent(IN)            :: file_id      ! File identifier
  integer, dimension(3), intent(IN)     :: field_dim    ! Vector of fiel dimension
  character(LEN=100), intent(IN)        :: field_name   ! Name of the data_set
  character(LEN=100), intent(IN)        :: unit         ! String identifying the unit
  real, intent(IN)                      :: conversion   ! Constant to multiply field with
  integer, intent(INOUT)                :: io
  real, dimension(field_dim(1),field_dim(2),field_dim(3)), intent(OUT) :: field ! Data to be stored

!-- local 

!--HDF5 Variables
  integer(HID_T) :: dset_id       ! Dataset identifier 
  integer(HID_T) :: filespace     ! Dataspace identifier in file 
  integer(HID_T) :: memspace      ! Dataspace identifier in memory
  integer(HID_T) :: plist_id      ! Property list identifier
  integer(HID_T) :: write_type    ! Type of data to write (REAL, DOUBLE, etc...)

! HDF5 variables for attribute
  character(LEN=100),parameter :: attribute_name1 = "UNITS" ! Attribute name
  character(LEN=100),parameter :: attribute_name2 = "CONV_COEFF" ! Attribute name

  integer(HID_T)  :: attr_id       ! Attribute identifier
  integer(HID_T)  :: aspace_id     ! Attribute Dataspace identifier
  integer(HID_T)  :: atype_id      ! Attribute Dataspace identifier
  integer(HSIZE_T), dimension(1) :: adims = (/1/) ! Attribute dimension
  integer         :: arank = 1                    ! Attribure rank
  integer(SIZE_T) :: attrlen       ! Length of the attribute string

  integer(HSIZE_T), allocatable :: dimsf(:) ! Dataset dimensions
                                            ! in the file.
  integer(HSIZE_T),  allocatable :: chunk_dims(:) ! Chunks dimensions
  integer(HSIZE_T),  allocatable :: count_ele(:)  
  integer(HSSIZE_T), allocatable :: offset(:) 
  integer(HSIZE_T) , allocatable :: stride(:)
  integer(HSIZE_T) , allocatable :: block(:)

  real, allocatable, target       :: reshape_field(:,:,:)
  real, pointer, dimension(:)     :: point_1D
  real, pointer, dimension(:,:)   :: point_2D
  real, pointer, dimension(:,:,:) :: point_3D

  integer :: rank
  integer :: i ! dummy index

  integer :: xid, yid, zid

  xid = myid_cart(1)
  yid = myid_cart(2)
  zid = myid_cart(3)

  ! Define rank of the array field

  rank = 0

  do i = 1, 3
    if (field_dim(i)>1) rank = rank + 1
  enddo

  if (rank == 0) then
    rank = 1
  endif

!----------------------------------------------------------------------------------------
  if (KIND(field) == 4) then      ! Get type of field, Could be real (Kind = 4) 
    write_type = H5T_NATIVE_REAL  ! or double (Kind = 8). Need for HDF5 subroutines
  else if(KIND(field) == 8) then
    write_type = H5T_NATIVE_DOUBLE
  else
    if (myid == 0) then
      write(io,*) 'The datatype to be read with HDF5 subroutine should be either REAL or DOUBLE.'
      write(io,*) 'ERROR in read_Field_dataset. Wrong datatype. Skip reading procedure.'         
      return
    endif
  endif

  ! Allocate variables needed for I/O
  allocate(dimsf(rank))
  allocate(chunk_dims(rank))
  allocate(count_ele(rank))
  allocate(offset(rank))
  allocate(stride(rank))
  allocate(block(rank))

  allocate(reshape_field(field_dim(3),field_dim(2),field_dim(1)))

  ! Define shape of dataset in file and memory
  ! Set up for 2D only, Need addition of extra case for 3D, i.e. case when 
  ! field_dims(3)/=1 when all the other are set to 1

  Select case(rank)
    case(3)

      dimsf      = (/nz_g, ny_g, nx_g/)
      chunk_dims = (/field_dim(3),field_dim(2),field_dim(1)/)
      point_3D   => reshape_field

      stride(1) = 1 
      stride(2) = 1
      stride(3) = 1 

      count_ele(1)  =  1 
      count_ele(2)  =  1
      count_ele(3)  =  1

      block(1)  = chunk_dims(1)
      block(2)  = chunk_dims(2)
      block(3)  = chunk_dims(3)

      offset(1) = zid*chunk_dims(1)
      offset(2) = yid*chunk_dims(2)
      offset(3) = xid*chunk_dims(3)
    case(2)  
      dimsf      = (/ny_g, nx_g/)
      chunk_dims = (/field_dim(2),field_dim(1)/)
      point_2D   => reshape_field(1,:,:)
      stride(1) = 1 
      stride(2) = 1 
  
      count_ele(1)  =  1 
      count_ele(2)  =  1 
      
      block(1)  = chunk_dims(1)
      block(2)  = chunk_dims(2)

      offset(1) = yid*chunk_dims(1)
      offset(2) = xid*chunk_dims(2)

    case(1)
      if (field_dim(1)/=1) then
        point_1D  => reshape_field(1,1,:)
        dimsf        = (/nx_g/)
        chunk_dims   = (/field_dim(1)/)
        count_ele(1) = 1
        offset(1)    = xid*chunk_dims(1)
        stride(1)    = 1
        block(1)     = chunk_dims(1)
      elseif (field_dim(2)/=1) then
        point_1D  => reshape_field(1,:,1)
        dimsf        = (/ny_g/)
        chunk_dims   = (/field_dim(2)/)
        count_ele(1) = 1
        offset(1)    = yid*chunk_dims(1)
        stride(1)    = 1
        block(1)     = chunk_dims(1)
      else
        point_1D  => reshape_field(:,1,1)
        dimsf        = (/nz_g/)
        chunk_dims   = (/field_dim(3)/)
        count_ele(1) = 1
        offset(1)    = zid*chunk_dims(1)
        stride(1)    = 1
        block(1)     = chunk_dims(1)
      endif
      
  end select

  ! allocate arrays with corresponding number of dimensions

  call h5screate_simple_f(rank, dimsf, filespace, hdf5_error)
  call h5screate_simple_f(rank, chunk_dims, memspace, hdf5_error)
  !
  ! Create chunked dataset.
  !
  call h5pcreate_f(H5P_DATASET_ACCESS_F, plist_id, hdf5_error)
!  call h5pset_chunk_f(plist_id, rank, chunk_dims, hdf5_error)

  call h5dopen_f(file_id, trim(field_name), dset_id, hdf5_error, plist_id)

!   ! Close data space created to free memory
!   !
   call h5sclose_f(filespace, hdf5_error)

  call h5dget_space_f(dset_id, filespace, hdf5_error)
  call h5sselect_hyperslab_f (filespace, H5S_SELECT_SET_F, offset, count_ele, hdf5_error, &
                                 stride, block)
  !
  ! Create property list for collective dataset read
  !
  call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdf5_error) 
  call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdf5_error)     
  !
  ! Write the dataset collectively. 
  !
  select case (rank)
    case(3)
      call h5dread_f(dset_id, write_type, point_3D, dimsf, hdf5_error, &
                    file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)
    case(2)
      call h5dread_f(dset_id, write_type, point_2D, dimsf, hdf5_error, &
                    file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)
    case(1)
      call h5dread_f(dset_id, write_type, point_1D, dimsf, hdf5_error, &
                    file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)
  end select

  ! Close dataspaces.
  call h5sclose_f(filespace, hdf5_error)
  call h5sclose_f(memspace, hdf5_error)
  ! Close the dataset.
  call h5dclose_f(dset_id, hdf5_error)
  ! Close the property list.
  call h5pclose_f(plist_id, hdf5_error)

  field = reshape(reshape_field,(/field_dim(1),field_dim(2),field_dim(3)/),ORDER=(/3,2,1/))       
  field = field/conversion

  select case (rank)
    case(3)
      nullify(point_3D)
    case(2)
      nullify(point_2D)
    case(1)
      nullify(point_1D)
   end select
  
  deallocate(dimsf)
  deallocate(chunk_dims)
  deallocate(count_ele)
  deallocate(offset)
  deallocate(stride)
  deallocate(block)
!----------------------------------------------------------------------------------------
  return
  end subroutine read_File_dataset

!=========================================================================================
  subroutine write_savefile(temp, pressure, div_qr,io)
!=========================================================================================
! writes data file for post-processing
! HDF5 format (one file for ALL the processes)
!-----------------------------------------------------------------------------------------

  use param_m,     only : x, y, z, nx, ny, nz, nsc, l_ref
  use topology_m,  only : myid, gcomm, ierr, MPI_INFO_NULL, hdf5_error

  use HDF5

  implicit none
!-----------------------------------------------------------------------------------------
! declarations passed in

  real, intent(in), dimension(:,:,:) :: temp
  real, intent(in), dimension(:,:,:) :: pressure
  real, intent(in), dimension(:,:,:) :: div_qr

  integer, intent(inout) :: io

! other declarations

! HDF5 PARAMETER DECLARATION
  integer(HID_T)     :: file_id       ! File identifier 
  integer(HID_T)     :: plist_id      ! Property list identifier
  integer(HID_T)     :: group_id      ! group identifier
  integer            :: error
  integer            :: info

  real, parameter    :: unity = 1.0

  character(LEN=100) :: field_name
  character(LEN=100) :: unit
  character(LEN=200) :: groupname
  integer, dimension(3) :: field_dim ! Vector of field dimension

  integer            :: i, j ! Dummy indeces

  character      :: getchar
  character*5    :: myid_ext
  character*7    :: i_time_ext
  character*11   :: time_ext
  character*100  :: filename, tarcmd, tartgt, tarsrc, dirname

  real           :: hr_conv

  logical :: exist


  if(myid.eq.0) then
    write(io,*) 'DTM: Writing save files.'
    write(io,*)
  endif

!----------------------------------------------------------------------------------------
! Use HDF5 file format to print the data

  ! Setup file access property list with parallel I/O access.
  info = MPI_INFO_NULL

  call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdf5_error)
  call h5pset_fapl_mpio_f(plist_id, gcomm, info, hdf5_error)
  !
  ! Create the file collectively.
  !
  filename = 'DTM_3D_MPI_SAVE.h5' 

  call h5fcreate_f(trim(filename), H5F_ACC_TRUNC_F, file_id, hdf5_error, &
                   access_prp = plist_id)

  call h5pclose_f(plist_id, hdf5_error)

!----------------------------------------------------------------------------------------
! Write the file attribute
  call write_run_attribute(file_id,io)

!----------------------------------------------------------------------------------------
! Write grid points

! Create group <grid>

  groupname = "Grid"
  call h5gcreate_f(file_id, groupname, group_id, hdf5_error)

  field_dim  = (/nx,1,1/)  
  field_name = 'x' 
  unit       = 'm'
  call write_File_dataset(group_id,x,field_dim,field_name,unit,l_ref, io)

  field_dim  = (/1,ny,1/)  
  field_name = 'y' 
  call write_File_dataset(group_id,y,field_dim,field_name,unit,l_ref, io)

  field_dim  = (/1,1,nz/)  
  field_name = 'z' 
  call write_File_dataset(group_id,z,field_dim,field_name,unit,l_ref, io)

  call h5gclose_f(group_id, hdf5_error)

  groupname = "Flow properties"
  call h5gcreate_f(file_id, groupname, group_id, hdf5_error)

  field_dim  = (/nx,ny,nz/)  
  field_name = 'Pressure' 
  unit       = 'Pa'
  call write_File_dataset(group_id,pressure,field_dim,field_name,unit,1.0, io)

  field_name = 'Temperature' 
  unit       = 'K'
  call write_File_dataset(group_id,temp,field_dim,field_name,unit,1.0, io)

  call h5gclose_f(group_id, hdf5_error)

  groupname = "Radiation Variables"
  call h5gcreate_f(file_id, groupname, group_id, hdf5_error)

!----------------------------------------------------------------------------------------
! Write variable Radiative source
  field_name = 'DivQr'
  unit       = ''
  hr_conv    = 1
  call write_File_dataset(group_id,div_qr,field_dim,field_name,unit,1.0, io)
  call h5gclose_f(group_id, hdf5_error)

!----------------------------------------------------------------------------------------
! Close HDF5 file
  call h5fclose_f(file_id, hdf5_error)

!-----------------------------------------------------------------------------------------
  return
  end subroutine write_savefile

!========================================================================================
  subroutine print_ray(i_ray, i_point, i_unit)
!========================================================================================
! Work only with 1 processors. 
! Print the segment from a given point in the file defined by i_unit
! FOR DEBUGGING PURPOSES ONLY
!----------------------------------------------------------------------------------------
  implicit none

  integer, intent(in) :: i_ray, i_point, i_unit

! local variables

  CHARACTER(LEN=256) :: name, Char_point, Char_ray
   
  write(Char_point,'(i6)') i_point
  write(Char_ray,'(i6)')   i_ray

  name = 'Point_' // trim(adjustl(Char_point)) // '_Ray_' // trim(adjustl(Char_ray))

  call tecplot_write_xy_line(i_unit, ray(i_ray,i_point),name)
!----------------------------------------------------------------------------------------
  return
  end subroutine print_ray

!========================================================================================
  subroutine tecplot_write_xy_line (iunit, ray_to_draw, name_ray)
!========================================================================================
!*******************************************************************************
!
!! TECPLOT_WRITE_XY_LINE writes out line data in XY geometry for use by TECPLOT.
!
!
!  Discussion:
!
!    The data is written as a GEOMETRY record.
!
!    The GEOMETRY record arguments, which you might want to adjust, include:
!
!      X = 0.0,
!      Y = 0.0, a fixed offset to be added to the (X,Y) data;
!      T = LINE, specifies that a line is being drawn;
!      CS = GRID, specifies that X and Y are measured in physical units;
!      L = SOLID, chooses a solid line type;
!      LT = 0.005, sets the line thickness, in frame units;
!      C = BLACK, chooses the line color to be black;
!      FC = BLACK, chooses the fill color, to be used if the line forms an area;
!      F = POINT, specifies that data will be (X1,Y1), (X2,Y2), ...(XN,YN);
!      S = GLOBAL, says that this line should appear in all "like" frames;
!
!    I am not completely clear on the filled areas.  In particular, I don't
!    understand whether:
!
!      A) the points are drawn, and if they enclose an area, the area is filled
!         with color FC;
!
!      B) the points are drawn, and the last point connected to the first,
!         and then, if FC is specified, that area is filled.
!
!      C) the points are drawn, and if FC is specified, the last point is
!         connected to the first and the area is filled.
!
!    The ZN parameter can be used to attact the line to a specific zone
!    or XY mapping.
!
!  Modified:
!
!    18 July 2001
!
!  Reference:
!
!    Section 3.5.2, "Geometry Record",
!    TECPLOT User's Manual, Version 7,
!    AMTEC Engineering, August 1996.
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IUNIT, the FORTRAN unit number associated with the file.
!
!    Input, integer NP, the number of nodes.
!
!    Input, real X(NP), the X coordinates of the nodes.
!
!    Input, real Y(NP), the Y coordinates of the nodes.
!
  implicit none

 ! Variable passed in:

  integer, intent(in)    :: iunit
  type(beam), intent(in) :: ray_to_draw

! local variables:

  integer :: n_path, i_path
  real    :: x, y, z
  CHARACTER(LEN=256) :: name_ray
 
  n_path = ray_to_draw%n_path


  write ( iunit, '(a)' ) 'TITLE = "Debug: Print Ray" '
  write ( iunit, '(a)' ) 'VARIABLES = "X (cm)", "Y (cm)", "Z (cm)","Temp_origin (K)"'
  write ( iunit, 10    )  trim(name_ray), n_path+2, n_path+1 

  do i_path = 1, n_path+1
    x = ray_to_draw%Point(i_path)%x
    y = ray_to_draw%Point(i_path)%y
    z = ray_to_draw%Point(i_path)%z
    write ( iunit, '(4g15.6)' ) x*1.0d+2, y*1.0d+2, z*1.0d+2, ray_to_draw%Temp_origin
  end do

  ! Print Origin point for debugging purposes

  x = ray_to_draw%Facet_origin%center%x
  y = ray_to_draw%Facet_origin%center%y
  z = ray_to_draw%Facet_origin%center%z
  write ( iunit, '(4g15.6)' ) x*1.0d+2, y*1.0d+2, z*1.0d+2,ray_to_draw%Facet_origin%Temp

  do i_path = 1, n_path + 1
    write ( iunit, '(2i5)' ) i_path, i_path + 1
  end do

!----------------------------------------------------------------------------------------
10 format ('ZONE T= "', (A), '" DATAPACKING=POINT, NODES=' , i5 ,', ELEMENTS=', i5 ,', ZONETYPE=FELINESEG')
  return
  end subroutine tecplot_write_xy_line

!========================================================================================
  subroutine Print_all_rays()
!========================================================================================
! Wrapper to print all rays. Use for debugging purposes

  ! Local variables

  integer :: i_ray, i_point, i_unit

! Create or Open file Ray.txt
  i_unit = 205
  open(i_unit, file='Ray.txt')

! Print all rays
  do i_point = 1, n_points
    do i_ray = 1, n_rays
      call print_ray(i_ray, i_point,i_unit)
    enddo
  enddo

! Close file 
  close(i_unit)  

!----------------------------------------------------------------------------------------
  return
  end subroutine Print_all_rays


!========================================================================================
  subroutine Test_rays()
!========================================================================================
! This subroutine tests all the rays constructed. Targets that the origin point is 
! consistent with the n_path^th point of the ray. if not, then the subroutine
! prints the defective ray. At the end of the tests, the subroutine prints some
! statistics. This function if for for debugging ONLY.
! Recall: 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top 
!----------------------------------------------------------------------------------------  

! Local variable
  integer :: i_unit, n_rays_tested
  integer :: i_dim1, i_dim2, i_dim3
  integer :: j_point, i_ray, i_boundary
  integer :: n_dim1, n_dim2, n_dim3

  integer, dimension(6) :: bad_rays_plane
  integer :: bad_rays_total
  integer :: n_rays_tested

  logical :: good_ray

  bad_rays_plane = 0
  bad_rays_total = 0 
  n_rays_tested  = 0

  i_unit = 205
  open(i_unit, file='Bad_Ray.txt')
  
  do i_boundary = 1, 6 

    select case(i_boundary)
      case(1:2)
        n_dim1 = mx-1
        n_dim2 = my-1
        n_dim3 = 1
      case(3:4)
        n_dim1 = 1
        n_dim2 = my-1
        n_dim3 = mz-1
      case(5:6)
        n_dim1 = mx-1
        n_dim2 = 1
        n_dim3 = mz-1
    end select
   
    do i_dim3 = 1, n_dim3
      do i_dim2 = 1, n_dim2
        do i_dim1 = 1, n_dim1

          j_point = Get_Points_Index(i_dim1,i_dim2,i_dim3,i_boundary)
          
          do i_ray = 1, n_rays

            n_rays_tested = n_rays_tested + 1
            good_ray      = consistent_ray(ray(i_ray,j_point))
                    
            if (.not.(good_ray)) then 
              bad_rays_total = bad_rays_total + 1 
              bad_rays_plane(i_boundary) = bad_rays_plane(i_boundary) + 1
              call print_ray(i_ray,j_point,i_unit)
            else
              cycle
            endif

          enddo
        enddo
      enddo
    enddo

  enddo

  close(i_unit)

! Write report

  write(*,*) ' '
  write(*,*) 'Subroutine Test_rays report'
  write(*,*) '==========================='
  write(*,*) n_rays_tested, ' rays tested out of ', n_rays*n_points
  write(*,*) 'There are ', bad_rays_total, ' bad rays, which is ', float(bad_rays_total)/float(n_rays_tested)*100.0, '%'
  write(*,*) 'There are ',  float(bad_rays_plane(1))/float(n_rays_tested)*100.0, ' % on South plane' 
  write(*,*) 'There are ',  float(bad_rays_plane(2))/float(n_rays_tested)*100.0, ' % on North plane'
  write(*,*) 'There are ',  float(bad_rays_plane(3))/float(n_rays_tested)*100.0, ' % on Right plane'
  write(*,*) 'There are ',  float(bad_rays_plane(4))/float(n_rays_tested)*100.0, ' % on Left plane'
  write(*,*) 'There are ',  float(bad_rays_plane(5))/float(n_rays_tested)*100.0, ' % on Bottom plane'
  write(*,*) 'There are ',  float(bad_rays_plane(6))/float(n_rays_tested)*100.0, ' % on Top plane'
  write(*,*) ' ' 
 
!----------------------------------------------------------------------------------------
  return
  end subroutine Test_rays 

!========================================================================================
  pure function Consistent_ray(Test_ray)
!========================================================================================
! Test if the origin point of the ray is consistent with the last point the ray cross.
! The last point the ray crosses should belong to the same plane as the origin.
! This function is for debugging purposes only as it requires a specific shape of the 
! type beam that keep in memory the coordinates of all the points the ray passes throught
!----------------------------------------------------------------------------------------

! Variables passed in

  type(beam), intent(in), target :: Test_ray

  logical :: Consistent_ray

! Local variables

  type(coordinates), pointer :: Point_origin, Last_point

! initialize logical Consistent_ray
  Consistent_ray = .false.

  if (associated(Last_point))   nullify(Last_point)
  if (associated(Point_origin)) nullify(Point_origin)
  
  Last_point   => Test_ray%Point(Test_ray%n_path+1)
  Point_origin => Test_ray%Facet_origin%center

  if (metric(Last_point%x,Point_origin%x) <= 1.0d-4) Consistent_ray = .true.
  if (metric(Last_point%y,Point_origin%y) <= 1.0d-4) Consistent_ray = .true.
  if (metric(Last_point%z,Point_origin%z) <= 1.0d-4) Consistent_ray = .true.

  if (associated(Last_point))   nullify(Last_point)
  if (associated(Point_origin)) nullify(Point_origin)
!----------------------------------------------------------------------------------------
  return
  end function Consistent_ray

!========================================================================================
  subroutine Print_irrad(i_proc)
!========================================================================================
! Print the vector irrad that contains the irradiation at the boundary condition
! Variables passed in:
! i_proc (integer,OPTIONAL) : index of the processor calling this subroutine
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in

  integer, intent(in), optional :: i_proc

! Local Variable
  real, parameter :: l_cgs = 100.0

  integer, target  :: i_dim1, i_dim2, i_dim3
  integer, pointer :: i_x, i_y

  integer :: n_x, n_y
  integer :: n_dim1, n_dim2, n_dim3

  integer :: i_unit, i_boundary
  integer :: j_point

  character(LEN=256) :: file_name

  i_unit = 1000

  if(present(i_proc)) then
    write(file_name,'(i5)') i_proc
    file_name = 'Irrad_' // trim(adjustl(file_name)) //'.dat'
  else
    file_name = 'Irrad.dat'
  endif

  open(i_unit, file=trim(file_name))
  write ( i_unit, '(a)' ) 'VARIABLES='
  write ( i_unit, '(a)' ) '"X (cm)"'
  write ( i_unit, '(a)' ) '"Y (cm)"'
  write ( i_unit, '(a)' ) '"Z (cm)"'
  write ( i_unit, '(a)' ) '"Irradiation (W/m2)"'
  write ( i_unit, '(a)' ) '"Radiosity (W/m2)"'
  write ( i_unit, '(a)' ) '"Temperature (K)"'

BOUNDARY: do i_boundary = 1, 6 

    if(associated(i_x)) nullify(i_x)
    if(associated(i_y)) nullify(i_y)

    select case(i_boundary)
      case(1:2)
        i_x => i_dim1
        i_y => i_dim2

        n_y = my-1
        n_x = mx-1

        n_dim1 = mx-1
        n_dim2 = my-1
        n_dim3 = 1
      case(3:4)
        i_x => i_dim2
        i_y => i_dim3

        n_y = mz-1
        n_x = my-1

        n_dim1 = 1
        n_dim2 = my-1
        n_dim3 = mz-1
      case(5:6)
        i_x => i_dim1
        i_y => i_dim3

        n_y = mz-1
        n_x = mx-1

        n_dim1 = mx-1
        n_dim2 = 1
        n_dim3 = mz-1
    end select
   
    select case(i_boundary)
      case(1)
          write ( i_unit,3) 'South', mx, my, 1
      case(2)
          write ( i_unit,3) 'North', mx, my, 1
      case(3)
          write ( i_unit,3) 'Right',  1, my, mz
      case(4)
          write ( i_unit,3) 'Left',   1, my, mz
      case(5)
          write ( i_unit,3) 'Top',   mx,  1, mz
      case(6)
          write ( i_unit,3) 'Bottom',mx,  1, mz
    end select

!----------------------------------------------------------------------------------------
! Write x

    do i_dim2 = 1, n_y, 1 
      do i_dim1 = 1, n_x, 1
        write (i_unit,'(1g15.6)') Boundaries(i_boundary)%element(i_dim1,i_dim2)%vertex(1)%x*l_cgs
      enddo
      write (i_unit,'(1g15.6)')   Boundaries(i_boundary)%element(n_x,i_dim2)%vertex(2)%x*l_cgs
    enddo 
  
    do i_dim1 = 1, n_x, 1
      write (i_unit,'(1g15.6)')   Boundaries(i_boundary)%element(i_dim1,n_y)%vertex(4)%x*l_cgs
    enddo
    write (i_unit,'(1g15.6)')     Boundaries(i_boundary)%element(n_x,n_y)%vertex(3)%x*l_cgs
 
!----------------------------------------------------------------------------------------
! Write y

    do i_dim2 = 1, n_y, 1 
      do i_dim1 = 1, n_x, 1
        write (i_unit,'(1g15.6)') Boundaries(i_boundary)%element(i_dim1,i_dim2)%vertex(1)%y*l_cgs
      enddo
      write (i_unit,'(1g15.6)')   Boundaries(i_boundary)%element(n_x,i_dim2)%vertex(2)%y*l_cgs
    enddo 
  
    do i_dim1 = 1, n_x, 1
      write (i_unit,'(1g15.6)')   Boundaries(i_boundary)%element(i_dim1,n_y)%vertex(4)%y*l_cgs
    enddo
    write (i_unit,'(1g15.6)')     Boundaries(i_boundary)%element(n_x,n_y)%vertex(3)%y*l_cgs


!----------------------------------------------------------------------------------------
! Write z

    do i_dim2 = 1, n_y, 1 
      do i_dim1 = 1, n_x, 1
        write (i_unit,'(1g15.6)') Boundaries(i_boundary)%element(i_dim1,i_dim2)%vertex(1)%z*l_cgs
      enddo
      write (i_unit,'(1g15.6)')   Boundaries(i_boundary)%element(n_x,i_dim2)%vertex(2)%z*l_cgs
    enddo 
  
    do i_dim1 = 1, n_x, 1
      write (i_unit,'(1g15.6)')   Boundaries(i_boundary)%element(i_dim1,n_y)%vertex(4)%z*l_cgs
    enddo
    write (i_unit,'(1g15.6)')     Boundaries(i_boundary)%element(n_x,n_y)%vertex(3)%z*l_cgs

!----------------------------------------------------------------------------------------
! Write Irradiation

    do i_dim3 = 1, n_dim3
      do i_dim2 = 1, n_dim2
        do i_dim1 = 1, n_dim1

          j_point = Get_Points_Index(i_dim1,i_dim2,i_dim3,i_boundary)
          write (i_unit,'(1g15.6)') irrad(j_point)

        enddo
      enddo
    enddo

!----------------------------------------------------------------------------------------
! Write Radiosity

    do i_dim3 = 1, n_dim3
      do i_dim2 = 1, n_dim2
        do i_dim1 = 1, n_dim1

          j_point = Get_Points_Index(i_dim1,i_dim2,i_dim3,i_boundary)
          write (i_unit,'(1g15.6)') pi*emission(j_point)

        enddo
      enddo
    enddo

!----------------------------------------------------------------------------------------
! Write Temperature

    do i_dim3 = 1, n_dim3
      do i_dim2 = 1, n_dim2
        do i_dim1 = 1, n_dim1
          write (i_unit,'(1g15.6)') Boundaries(i_boundary)%element(i_x,i_y)%Temp
        enddo
      enddo
    enddo

    if(associated(i_x)) nullify(i_x)
    if(associated(i_y)) nullify(i_y)

  enddo BOUNDARY

  close(i_unit)

  return

3  format(' ZONE t="',A8,'", i=',i5,', j=',i5,', k=',i5,', DATAPACKING=block, VARLOCATION=([4,5,6]=CELLCENTERED) DT=(DOUBLE,DOUBLE,DOUBLE, DOUBLE, DOUBLE, DOUBLE)')
!----------------------------------------------------------------------------------------
  end subroutine Print_irrad

!========================================================================================
  subroutine Print_Fluxes(myid,io)
!========================================================================================
! Computes the different integrated radiative fluxes (Incident, leaving, net)
! and print them for each faces

! Variables passed in
!--------------------
! myid: id of the processor calling the subroutine
! io: integer, unit of the output file provided by the user
!----------------------------------------------------------------------------------------

  implicit none

  include 'mpif.h'

! Variables passed in 

  integer, intent(in) :: myid
  integer, intent(in) :: io

! Local variables

  real, allocatable, dimension(:,:) :: Array1, Array2

  real, dimension(6) :: Incident,   Leaving,   Net_flux
  real, dimension(6) :: Incident_g, Leaving_g, Net_flux_g ! (Global arrays, to be printed by process 0)

  real, parameter :: l_ref = 1.0

  integer :: jx, jy, jz
  integer :: i_dimension

  integer, dimension(1:2) :: order1 = (/2,1/)

  integer :: lb, ub

!------------------------------------------------------------------------------
! South

  allocate(Array1(mx-1,my-1))
  allocate(Array2(mx-1,my-1))
  Array1 = 0.0 
  Array2 = 0.0

! Use reshape array

  lb = Get_Points_Index(1,1,1,1)
  ub = Get_Points_Index(mx-1,my-1,1,1) 

  Array1 = RESHAPE(irrad(lb:ub),(/mx-1,my-1/),ORDER = ORDER1)
  Array2 = RESHAPE(pi*emission(lb:ub),(/mx-1,my-1/),ORDER = ORDER1)

  Incident(1) = 0.0
  Leaving(1)  = 0.0
  Net_flux(1) = 0.0

  do jy = 1, my-1
    do jx = 1, mx-1
      Incident(1) = Incident(1) + Boundaries(1)%element(jx,jy)%Surface*Array1(jx,jy)
      Leaving(1)  = Leaving(1)  + Boundaries(1)%element(jx,jy)%Surface*Array2(jx,jy)
    enddo
  enddo 

  Incident(1) = Incident(1)*l_ref**2
  Leaving(1)  = Leaving(1)*l_ref**2
  Net_flux(1) = Leaving(1) - Incident(1)

  deallocate(Array1)
  deallocate(Array2)

!------------------------------------------------------------------------------
! North

  allocate(Array1(mx-1,my-1))
  allocate(Array2(mx-1,my-1))
  Array1 = 0.0 
  Array2 = 0.0

! Use reshape array

  lb = Get_Points_Index(1,1,mz,2)
  ub = Get_Points_Index(mx-1,my-1,mz,2) 

  Array1 = RESHAPE(irrad(lb:ub),(/mx-1,my-1/),ORDER = ORDER1)
  Array2 = RESHAPE(pi*emission(lb:ub),(/mx-1,my-1/),ORDER = ORDER1)

  Incident(2) = 0.0
  Leaving(2)  = 0.0
  Net_flux(2) = 0.0

  do jy = 1, my-1
    do jx = 1, mx-1
      Incident(2) = Incident(2) + Boundaries(2)%element(jx,jy)%Surface*Array1(jx,jy)
      Leaving(2)  = Leaving(2)  + Boundaries(2)%element(jx,jy)%Surface*Array2(jx,jy)
    enddo
  enddo 

  Incident(2) = Incident(2)*l_ref**2
  Leaving(2)  = Leaving(2)*l_ref**2
  Net_flux(2) = Leaving(2) - Incident(2)

  deallocate(Array1)
  deallocate(Array2)

!------------------------------------------------------------------------------
! Right plane: (y,z,x)

  allocate(Array1(my-1,mz-1))
  allocate(Array2(my-1,mz-1))
  Array1 = 0.0 
  Array2 = 0.0

! Use reshape array
  lb = Get_Points_Index(1,1,1,3)
  ub = Get_Points_Index(1,my-1,mz-1,3) 

  Array1 = RESHAPE(irrad(lb:ub),(/my-1,mz-1/),ORDER = ORDER1)
  Array2 = RESHAPE(pi*emission(lb:ub),(/my-1,mz-1/),ORDER = ORDER1)

  Incident(3) = 0.0
  Leaving(3)  = 0.0
  Net_flux(3) = 0.0

  do jz = 1, mz-1
    do jy = 1, my-1
      Incident(3) = Incident(3) + Boundaries(3)%element(jy,jz)%Surface*Array1(jy,jz)
      Leaving(3)  = Leaving(3)  + Boundaries(3)%element(jy,jz)%Surface*Array2(jy,jz)
    enddo
  enddo 

  Incident(3) = Incident(3)*l_ref**2
  Leaving(3)  = Leaving(3)*l_ref**2
  Net_flux(3) = Leaving(3) - Incident(3)

  deallocate(Array1)
  deallocate(Array2)

!------------------------------------------------------------------------------
! Left plane: (-y,z,-x)

  allocate(Array1(my-1,mz-1))
  allocate(Array2(my-1,mz-1))
  Array1 = 0.0 
  Array2 = 0.0

! Use reshape array
  lb = Get_Points_Index(mx,   1,   1, 4)
  ub = Get_Points_Index(mx,my-1,mz-1, 4) 

  Array1 = RESHAPE(irrad(lb:ub),(/my-1,mz-1/),ORDER = ORDER1)
  Array2 = RESHAPE(pi*emission(lb:ub),(/my-1,mz-1/),ORDER = ORDER1)

  Incident(4) = 0.0
  Leaving(4)  = 0.0
  Net_flux(4) = 0.0

  do jz = 1, mz-1
    do jy = 1, my-1
      Incident(4) = Incident(4) + Boundaries(4)%element(jy,jz)%Surface*Array1(jy,jz)
      Leaving(4)  = Leaving(4)  + Boundaries(4)%element(jy,jz)%Surface*Array2(jy,jz)
    enddo
  enddo 

  Incident(4) = Incident(4)*l_ref**2
  Leaving(4)  = Leaving(4)*l_ref**2
  Net_flux(4) = Leaving(4) - Incident(4)

  deallocate(Array1)
  deallocate(Array2)

!------------------------------------------------------------------------------
! Bottom plane: (z,x,y)

  allocate(Array1(mz-1,mx-1))
  allocate(Array2(mz-1,mx-1))
  Array1 = 0.0 
  Array2 = 0.0

! Use reshape array
  lb = Get_Points_Index(1   ,   1,    1, 5)
  ub = Get_Points_Index(mx-1,   1, mz-1, 5) 

  Array1 = RESHAPE(irrad(lb:ub),(/mz-1,mx-1/),ORDER = ORDER1)
  Array2 = RESHAPE(pi*emission(lb:ub),(/mz-1,mx-1/),ORDER = ORDER1)

  Incident(5) = 0.0
  Leaving(5)  = 0.0
  Net_flux(5) = 0.0

  do jz = 1, mz-1
    do jx = 1, mx-1
      Incident(5) = Incident(5) + Boundaries(5)%element(jx,jz)%Surface*Array1(jz,jx)
      Leaving(5)  = Leaving(5)  + Boundaries(5)%element(jx,jz)%Surface*Array2(jz,jx)
    enddo
  enddo 

  Incident(5) = Incident(5)*l_ref**2
  Leaving(5)  = Leaving(5)*l_ref**2
  Net_flux(5) = Leaving(5) - Incident(5)

  deallocate(Array1)
  deallocate(Array2)

!------------------------------------------------------------------------------
! Top plane (z,-x,-y)

  allocate(Array1(mz-1,mx-1))
  allocate(Array2(mz-1,mx-1))
  Array1 = 0.0 
  Array2 = 0.0

! Use reshape array
  lb = Get_Points_Index(1   ,   1,    1, 6)
  ub = Get_Points_Index(mx-1,   1, mz-1, 6) 

  Array1 = RESHAPE(irrad(lb:ub),(/mz-1,mx-1/),ORDER = ORDER1)
  Array2 = RESHAPE(pi*emission(lb:ub),(/mz-1,mx-1/),ORDER = ORDER1)

  Incident(6) = 0.0
  Leaving(6)  = 0.0
  Net_flux(6) = 0.0

  do jz = 1, mz-1
    do jx = 1, mx-1
      Incident(6) = Incident(6) + Boundaries(6)%element(jx,jz)%Surface*Array1(jz,jx)
      Leaving(6)  = Leaving(6)  + Boundaries(6)%element(jx,jz)%Surface*Array2(jz,jx)
    enddo
  enddo 

  Incident(6) = Incident(6)*l_ref**2
  Leaving(6)  = Leaving(6)*l_ref**2
  Net_flux(6) = Leaving(6) - Incident(6)

  deallocate(Array1)
  deallocate(Array2)

!----------------------------------------------------------------------------------------
! If the face is a virtual boundary, then set the Incident fluxes, Leaving flux, and 
! net_flux to 0

  do i_dimension = 1, 6
    if (is_vb(i_dimension)) then 
      Incident(i_dimension) = 0.0d+0
      Leaving(i_dimension)  = 0.0d+0
      Net_flux(i_dimension) = 0.0d+0
    endif
  enddo

!----------------------------------------------------------------------------------------
! For each face, get the sum from all the processors. Print the result on the master
! process (process id 0)

  do i_dimension = 1, 6

! Incident flux
    call MPI_REDUCE(Incident(i_dimension),Incident_g(i_dimension),1,MPI_DOUBLE_PRECISION,&
                    MPI_SUM,0,cart_comm,ierr_MPI_rad)
! Radiosity flux
    call MPI_REDUCE(Leaving(i_dimension),Leaving_g(i_dimension)  ,1,MPI_DOUBLE_PRECISION,&
                    MPI_SUM,0,cart_comm,ierr_MPI_rad)
! Net flux
    call MPI_REDUCE(Net_flux(i_dimension),Net_flux_g(i_dimension),1,MPI_DOUBLE_PRECISION,&
                    MPI_SUM,0,cart_comm,ierr_MPI_rad)

  enddo

!----------------------------------------------------------------------------------------
! Master process (myid==0) prints the results in file given by unit io

  if (myid==0) then
    write(io,*) char(10)
    write(io,*) '========================================================================='
    write(io,*) 'Print Integrated Fluxes over each face of the domain'
    write(io,*) '========================================================================='
    write(io,*) 'FACE:     Irradiation (W),    Radiosity (W),     Net_flux (W)'
    write(io,*) 'South  ', char(9), Incident_g(1), char(9), Leaving_g(1), char(9), Net_flux_g(1)
    write(io,*) 'North  ', char(9), Incident_g(2), char(9), Leaving_g(2), char(9), Net_flux_g(2)
    write(io,*) 'Right  ', char(9), Incident_g(3), char(9), Leaving_g(3), char(9), Net_flux_g(3)
    write(io,*) 'Left   ', char(9), Incident_g(4), char(9), Leaving_g(4), char(9), Net_flux_g(4)
    write(io,*) 'Bottom ', char(9), Incident_g(5), char(9), Leaving_g(5), char(9), Net_flux_g(5)
    write(io,*) 'Top    ', char(9), Incident_g(6), char(9), Leaving_g(6), char(9), Net_flux_g(6)
    write(io,*) '---------------------------------------------------------'
    write(io,*) 'Sum Net Flux (W):', (char(9),jx=1,7), sum(Net_flux_g)
    write(io,*) char(10)
  endif

!----------------------------------------------------------------------------------------
  return
  end subroutine Print_Fluxes
  
!========================================================================================
  subroutine Print_Div_qr(Div_qr)
!========================================================================================
! Print Divergence of qr in Tecplot file
!----------------------------------------------------------------------------------------

  implicit none

  real, intent(in) :: Div_qr(:,:,:)

  integer jx, jy, jz
  integer i_unit

  real, parameter :: l_cgs = 100.0

  i_unit = 1002

  open(i_unit, file='Div_qr.dat')
  write ( i_unit, '(a)' ) 'VARIABLES='
  write ( i_unit, '(a)' ) '"X (cm)"'
  write ( i_unit, '(a)' ) '"Y (cm)"'
  write ( i_unit, '(a)' ) '"Z (cm)"'
  write ( i_unit, '(a)' ) '"Div_qr (W/m3)"'
  write ( i_unit,3) 'Inside', mx, my, mz

  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') xx(jx)*l_cgs
      enddo
    enddo 
  enddo

  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') yy(jy)*l_cgs
      enddo
    enddo 
  enddo
 
  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') zz(jz)*l_cgs
      enddo
    enddo 
  enddo 

  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') Div_qr(jx,jy,jz)
      enddo
    enddo 
  enddo 

  close(i_unit)

3  format(' ZONE t="',A8,'", i=',i5,', j=',i5,', k=',i5,', DATAPACKING=block, DT=(DOUBLE,DOUBLE,DOUBLE,DOUBLE)')
!----------------------------------------------------------------------------------------
  return
  end subroutine Print_Div_qr

!========================================================================================
  subroutine Print_Cell_Source(i_proc)
!========================================================================================
! Print Divergence of qr in Tecplot file
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in

  integer, intent(in), optional :: i_proc

! Local Variables

  integer :: jx, jy, jz
  integer :: i_unit

  real, parameter :: l_cgs = 100.0

  character(LEN=256) :: file_name

  i_unit = 1003

  if(present(i_proc)) then
    write(file_name,'(i5)') i_proc
    file_name = 'Cell_source_' // trim(adjustl(file_name)) //'.dat'
  else
    file_name = 'Cell_source.dat'
  endif

  open(i_unit, file=trim(file_name))

  write ( i_unit, '(a)' ) 'VARIABLES='
  write ( i_unit, '(a)' ) '"X (cm)"'
  write ( i_unit, '(a)' ) '"Y (cm)"'
  write ( i_unit, '(a)' ) '"Z (cm)"'
  write ( i_unit, '(a)' ) '"Cell source (W)"'
  write ( i_unit, '(a)' ) '"Cell Temp (K)"'
  write ( i_unit, '(a)' ) '"Cell Abs (m-1)"'
  write ( i_unit, '(a)' ) '"Cell Surface (m2)"'
  write ( i_unit,3) 'Inside', mx, my, mz

  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') xx(jx)*l_cgs
      enddo
    enddo 
  enddo

  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') yy(jy)*l_cgs
      enddo
    enddo 
  enddo
 
  do jz = 1, mz
    do jy = 1, my
      do jx = 1, mx
	write (i_unit,'(1g15.6)') zz(jz)*l_cgs
      enddo
    enddo 
  enddo 

  do jz = 1, mz-1
    do jy = 1, my-1
      do jx = 1, mx-1
	write (i_unit,'(1g15.6)') cell_tot(jx,jy,jz)%Source_cv
      enddo
    enddo 
  enddo 

  do jz = 1, mz-1
    do jy = 1, my-1
      do jx = 1, mx-1
	write (i_unit,'(1g15.6)') cell_tot(jx,jy,jz)%Temp
      enddo
    enddo 
  enddo 

  do jz = 1, mz-1
    do jy = 1, my-1
      do jx = 1, mx-1
	write (i_unit,'(1g15.6)') cell_tot(jx,jy,jz)%c_abs
      enddo
    enddo 
  enddo 

  do jz = 1, mz-1
    do jy = 1, my-1
      do jx = 1, mx-1
	write (i_unit,'(1g15.6)') cell_tot(jx,jy,jz)%Surface
      enddo
    enddo 
  enddo 

  close(i_unit)

3  format(' ZONE t="',A8,'", i=',i5,', j=',i5,', k=',i5,', DATAPACKING=block,  VARLOCATION=([4,5,6,7]=CELLCENTERED) DT=(DOUBLE,DOUBLE,DOUBLE,DOUBLE)')
!----------------------------------------------------------------------------------------
  return
  end subroutine Print_Cell_Source

!========================================================================================
  subroutine Print_Diag_North
!========================================================================================
! Print the vector irrad that contains the irradiation at the boundary condition
!----------------------------------------------------------------------------------------

  implicit none

  real, parameter :: l_cgs = 100.0

  real, allocatable :: Array_to_print(:,:)
  real, allocatable :: Array2(:,:)
  real emitted_flux

  integer jx, jy
  integer i_unit
  integer, dimension(1:2) :: order1 = (/2,1/)

  integer lb, ub

  i_unit = 1000

  open(i_unit, file='Configuration_Factor.dat')
  write ( i_unit, '(a)' ) 'VARIABLES='
  write ( i_unit, '(a)' ) '"X (m)"'
  write ( i_unit, '(a)' ) '"Y (m)"'
  write ( i_unit, '(a)' ) '"Z (m)"'
  write ( i_unit, '(a)' ) '"View Factor"'
  write ( i_unit,3) 'North', mx-1, 1, 1

  do jx = 1, mx-1
    write (i_unit,'(1g15.6)') 0.5*(cell_tot(jx,jx,mz-1)%x(2)+cell_tot(jx,jx,mz-1)%x(1))
  enddo

  do jx = 1, mx-1
    write (i_unit,'(1g15.6)') 0.5*(cell_tot(jx,jx,mz-1)%y(2)+cell_tot(jx,jx,mz-1)%y(1))
  enddo

  do jx = 1, mx-1
    write (i_unit,'(1g15.6)') zz(mz)
  enddo

  allocate(Array_to_print(mx-1,my-1))
  allocate(Array2(mx-1,my-1))

  Array_to_print = 0.0 
  Array2         = 0.0

! Use reshape array

  lb = (mx-1)*(my-1)+1
  ub = 2*(mx-1)*(my-1)

  Array_to_print = RESHAPE(irrad(lb:ub),(/mx-1,my-1/),ORDER = ORDER1)

  lb = 1 
  ub = (mx-1)*(my-1)

  Array2 = RESHAPE(pi*emission(lb:ub),(/mz-1,mx-1/),ORDER = ORDER1)

! Compute the total flux emitted from the south face

  emitted_flux = 0.0

  do jy = 1, my-1
    do jx = 1, mx-1
      emitted_flux = emitted_flux + abs((xx(jx+1)-xx(jx))*(yy(jy+1)-yy(jy)))*Array2(jx,jy)
    enddo
  enddo 

  do jx = 1, mx-1
    write (i_unit,'(1g15.6)') Array_to_print(jx,jx)/emitted_flux
  enddo


! View factor: by definition, the view factor between two infinitesimal surface elements
! dAi and dAj is defined as:
!
! dF(Ai)-(Aj) = (diffuse energy leaving dAi directly toward and intercepted by dAj)/(total diffuse energy leaving dAi)

  deallocate(Array_to_print)
!------------------------------------------------------------------------------

  close(i_unit)

  return

3  format(' ZONE t="',A8,'", i=',i5,', j=',i5,', k=',i5,', DATAPACKING=block, DT=(DOUBLE,DOUBLE,DOUBLE, DOUBLE, DOUBLE)')
!----------------------------------------------------------------------------------------
  end subroutine Print_Diag_North

!========================================================================================
  subroutine Error_MSG(arg)
!========================================================================================
! Thiis program is private. Prints message error depending on the value of 
! arg.
! 
! Variables passed in:
!---------------------
! arg : integer. Indicate nature of error

! The unit for write operation is given by global variable output_unit
!
!----------------------------------------------------------------------------------------

  implicit none

! Variables passed in

  integer :: arg

! Local variables

  character(80) :: MSG

  select case (arg)
    case(1)
      MSG = 'ERROR!! Shapes of the first and second variable do not conform.'
    case(2)
      MSG = 'ERROR!! Shapes of the first and third variable do not conform.'
    case(3)
      MSG = 'ERROR in Get_facet_index. Boundary not assigned to facet.'
    case(4)
      MSG = 'ERROR in Get_facet_index. i_fix_coord of Boundary unknown.'
    case(5)
      MSG = 'ERROR in Get_Cell_Index. The point is not inside the domain.'
  end select

  write(output_unit,*) trim(MSG)  

  stop 
!----------------------------------------------------------------------------------------
  end subroutine Error_MSG
!----------------------------------------------------------------------------------------
  end module radiation_m


!========================================================================================
  module dtm_m
!========================================================================================
  implicit none

!  private

  type ribv_pointer
    real, pointer, dimension(:,:,:) :: array
  end type ribv_pointer

!------------------Define some parameters------------------------------------------------
  integer, parameter :: n_iteration = 5000   ! Maximum number of iteration for 
                                             ! RTE integration
  real,    parameter :: tolerance   = 1.0e-5 !tolerance for error_g 
                                             !(decrease it to achieve better solution)
!-----------------Define some type-------------------------------------------------------

! PUBLIC Procedures
!  public :: DTM 

! PUBLIC Variables

  contains

!========================================================================================
  subroutine DTM(temp,pressure,yspecies,div_qr,io)
!========================================================================================
! This version is for validation tests only.
! MPI Integrated subroutine
! Works for non-uniform grid 
!----------------------------------------------------------------------------------------
! First version (2D): 12/18/03
! by Yi Wang (yiwang@wam.umd.edu)
!----------------------------------------------------------------------------------------
! Adapted by Vivien Lecoustre (vlecous1@umd.edu)
! Successful tests periodic boundary conditions
! 07/25/12 
!=========================================================================================
! Variables passed IN
! temp     :: Real. Array of local temperature.     UNIT = Kelvin. 
! pressure :: Real. Array of local pressure.        UNIT = Pascals
! Yspecie  :: Real. Array of local species mass fraction. UNIT = Dimensionless
! io       :: Integer, optional. Refers to the units of the information output file.
!
! Variables passed OUT
! div_qr :: Real. Array containing the divergence of the local radiative heat flux. 
! UNIT = Dimensionless
! 
! if periodic_x = 1, then RTE has periodic BC in x direction
! if periodic_y = 1, then RTE has periodic BC in y direction 
! if periodic_z = 1, then RTE has periodic BC in z direction

!  use param_m, only : nx, ny, nz, nsc, periodic_x, periodic_y, periodic_z

  use radiation_m
  use topology_m, only : myid
  implicit none

  include 'mpif.h'

! Declarations passed in
  
  real, intent(in), dimension(:,:,:)   :: temp
  real, intent(in), dimension(:,:,:)   :: pressure
  real, intent(in), dimension(:,:,:,:) :: yspecies

  integer, intent(in), optional :: io                       ! Unit of the error output file                     

! Variables passed out
  real, intent(out), dimension(:,:,:), allocatable :: div_qr

! Local declarations

  type(ribv_pointer), dimension(6) :: New_rivb
  type(ribv_pointer), dimension(6) :: Old_rivb

  integer :: i_plane, i_point, i_iter
 
  integer :: nx, ny, nz

  integer, dimension(3) :: dim_xyz

  integer, dimension(2,6), target :: indices_plane ! vector containing the start and 
                                                  ! end indices for the planes

  logical, dimension(6) :: vb_flag      ! TRUE = virtual boundary, FALSE = real boundary
  
  integer, pointer :: i_1, i_2                     ! Pointers to be used with indices_plane

  real :: max_variation
  real :: error, error_g                      
  
!----------------------------------------------------------------------------------------
! If div_qr is allocated deallocate and reallocate

  if (allocated(div_qr)) deallocate(div_qr)

!----------------------------------------------------------------------------------------
! Ensure the variable have same dimensions

  if (all(shape(temp)/=shape(pressure))) then
    if (present(io)) then 
      call Error_MSG(1,io)
    else
      call Error_MSG(1)
    endif

  elseif (all(shape(temp)/=shape(yspecies(:,:,:,1)))) then
    if (present(io)) then 
      call Error_MSG(2,io)
    else
      call Error_MSG(2)
    endif

  else
    dim_xyz = shape(temp)
    nx = dim_xyz(1)
    ny = dim_xyz(2)
    nz = dim_xyz(3)
    allocate(div_qr(nx,ny,nz))
  endif

!----------------------------------------------------------------------------------------
! Create indices vector

! SouthernPoints
  indices_plane(1,1) = 1
  indices_plane(2,1) = (mx-1)*(my-1)
! NorthernPoints:
  indices_plane(1,2) = indices_plane(2,1)+1
  indices_plane(2,2) = indices_plane(2,1)+(mx-1)*(my-1)
! RightPoints
  indices_plane(1,3) = indices_plane(2,2)+1
  indices_plane(2,3) = indices_plane(2,2)+(my-1)*(mz-1)
! LeftPoints: 
  indices_plane(1,4) = indices_plane(2,3)+1
  indices_plane(2,4) = indices_plane(2,3)+(my-1)*(mz-1)
! BottomPoints: 
  indices_plane(1,5) = indices_plane(2,4)+1 
  indices_plane(2,5) = indices_plane(2,4)+(mx-1)*(mz-1)
! TopPoints
  indices_plane(1,6) = indices_plane(2,5)+1 
  indices_plane(2,6) = indices_plane(2,5)+(mx-1)*(mz-1)

!----------------------------------------------------------------------------------------
! Create vb_flag vector

  vb_flag = is_vb

!----------------------------------------------------------------------------------------
! Assign the pointers Old_rivb and New_rivb to their corresponding plan

! SouthernPoints
  Old_rivb(1)%array => rivb0_s
  New_rivb(1)%array => rivb_s
! NorthernPoints
  Old_rivb(2)%array => rivb0_n
  New_rivb(2)%array => rivb_n
! RightPoints
  Old_rivb(3)%array => rivb0_r
  New_rivb(3)%array => rivb_r
! LeftPoints
  Old_rivb(4)%array => rivb0_l
  New_rivb(4)%array => rivb_l
! BottomPoints
  Old_rivb(5)%array => rivb0_b
  New_rivb(5)%array => rivb_b
! TopPoints
  Old_rivb(6)%array => rivb0_t
  New_rivb(6)%array => rivb_t

! Refresh values of rivb and rivb0_b

  do i_plane = 1, 6
    Old_rivb(i_plane)%array = 0.0
    New_rivb(i_plane)%array = 0.0
  enddo

!----------------------------------------------------------------------------------------
! Initialization of div_qr values
  div_qr = 0.0

!----------------------------------------------------------------------------------------
! Update control volume properties: Temperature, Kappa, BlackBody emission

  call Update_Cell(yspecies,pressure,temp)

!----------------------------------------------------------------------------------------
! integrate the RTE along each ray
!  Recall: 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top
!---------------------------------------------------------------------------------------- 
 
  Iteration: do i_iter = 1, n_iteration    ! begin iteration

 ! Set variables to zero between each iteration
    error              = 0.0
    error_g            = 0.0
    cell_tot%Source_cv = 0.0
    irrad              = 0.0

!----------------------------------------------------------------------------------------
! Compute Radiative emission from boundary condition
    call Radiative_emission(temp,irrad0,emission)

!----------------------------------------------------------------------------------------
! Get virtual boundary Irradiation
    call exchange_RIVB   !rivb doesn't change from the last time integration

!----------------------------------------------------------------------------------------
! Loop over the boundaries
! Recall: 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top

    BOUNDARY: do i_plane = 1,6,1

      i_1 => indices_plane(1,i_plane)
      i_2 => indices_plane(2,i_plane)

      do i_point = i_1,i_2
        call int_RTE(i_point,i_plane,vb_flag(i_plane))
      enddo

      if (.not.(vb_flag(i_plane))) then 
        call get_error_1D(irrad0(i_1:i_2),irrad(i_1:i_2),max_variation,io)
      else
        call get_error_3D(Old_rivb(i_plane)%array,New_rivb(i_plane)%array,max_variation,io)
      endif 
    
      error = max(error,max_variation) 

    enddo BOUNDARY

!----------------------------------------------------------------------------------------
! Get the maximum value of error across the MPI communicator

    call MPI_ALLREDUCE(error,error_g,1,MPI_REAL8,MPI_MAX,cart_comm,ierr_MPI_rad)

    if (myid==0) write(io,*) 'Max_err', error_g
    if (error_g.lt.tolerance) exit

  enddo iteration

! print number of iteration
  if (myid==0) write(io,*) 'number of iteration for DTM is', i_iter 

!----------------------------------------------------------------------------------------
! Iteration finished
! Compute the divergence of radiative source terms

  call Get_Div_qr(div_qr)

! Recall: div_qr is in units of [W/m^3]. Length was computed from product [x,y,z]*l_ref

! Nondimentionalization of the source term div_qr
!
! div_qr=div_qr/(rho_ref*a_ref**3/l_ref)
!
!----------------------------------------------------------------------------------------
! Filter source term 

! call filter(div_qr,io)

!----------------------------------------------------------------------------------------
! Deallocate and nullify pointer
  do i_plane = 1, 6
    nullify(Old_rivb(i_plane)%array)
    nullify(New_rivb(i_plane)%array)
  enddo

  return
!------------------------------------------------------------------------------------------
  end subroutine DTM

!========================================================================================
  subroutine int_RTE(j_point, edge, vb_flag)
!========================================================================================
! integrate the RTE for all the rays incident on one bounday points
!----------------------------------------------------------------------------------------

  use radiation_m
 
  implicit none

! Declarations passed in

  integer,intent(in):: j_point   ! Index of number of boundary points 
  integer,intent(in):: edge      ! 1-South, 2-North, 3-Right, 4-Left, 5-Bottom, 6-Top 
  logical,intent(in):: vb_flag   ! True: virtual boundary, False: real boundary

! Local declarations

  real    :: tau_dir ! transmissivity along ray crossing a cuboid. It is introduced to
                     ! save CPU time

  integer :: m_x, m_y, m_z
  integer :: i_path, n_path_ray
  integer :: j_point_wall
  integer :: i_phi, i_theta
  integer :: i_ray, i_ray2
  integer :: target_origin

  integer, dimension(4) :: Index_Facet

  real :: inten_ghost
  real :: inten0, inten1, inten_b
  real :: surf_area
  
  logical :: vb_origin ! If vb_origin ==.TRUE. then the origin of the ray is on a virtual
                       ! boundary. Expression for the Irradiance must me changed

  real, pointer :: Array_incident(:,:,:)
  type(beam),   pointer :: Point_origin
  type(cuboid), pointer :: Cell_path

!----------------------------------------------------------------------------------------
! Associate pointer 'Array_incident' with target: Face of rib
! This is based on the value of integer edge.
! Recall: edge = 1 => South
!         edge = 2 => North
!         edge = 3 => Right
!         edge = 4 => Left
!         edge = 5 => Bottom
!         edge = 6 => Top

  Select Case(edge)
    case(1)  ! => South
      Array_incident => rib%South
    case(2)  ! => North
      Array_incident => rib%North
    case(3)  ! => Right
      Array_incident => rib%Right
    case(4)  ! => Left
      Array_incident => rib%Left
    case(5)  ! => Bottom
      Array_incident => rib%Bottom
    case(6)  ! => Top
      Array_incident => rib%Top
    case default
      Print *,' Error Variable edge cannoy be greater than 6'
      Print *,' Terminate execution'
      call allocate_radiation_arrays(-1)
      stop
  end select

!----------------------------------------------------------------------------------------
! Loop over rays

  Rays: do i_ray = 1, n_rays

    inten_ghost = 0.0

! Assign pointer Point_origin to origin(i_ray,j_point) for ease of manipulation

    Point_origin   => ray(i_ray,j_point) 

! Use Get_Facet_Coordinates to populate m_x, m_y, m_z, target_origin
    Index_Facet    = Get_Facet_Coordinates(Point_origin%Facet_origin) 

    m_x            = Index_Facet(1)
    m_y            = Index_Facet(2)
    m_z            = Index_Facet(3)
    target_origin  = Index_Facet(4)

    i_ray2         = Point_origin%iray2
    i_phi          = Point_origin%iphi
    i_theta        = Point_origin%itheta
    n_path_ray     = Point_origin%n_path

! If virtual voundaries or periodic conditions in x, y, z

    if ((target_origin==1).and.is_vb(1)) then 
      inten_ghost = rivb_s(i_ray2,m_x,m_y)  ! South plan
      vb_origin = .TRUE.
    elseif ((target_origin==2).and.is_vb(2)) then 
      inten_ghost = rivb_n(i_ray2,m_x,m_y)  ! North plan
      vb_origin = .TRUE.
    elseif ((target_origin==3).and.is_vb(3)) then 
      inten_ghost = rivb_r(i_ray2,m_y,m_z)  ! Right plan
      vb_origin = .TRUE.
    elseif ((target_origin==4).and.is_vb(4)) then 
      inten_ghost = rivb_l(i_ray2,m_y,m_z)  ! Left plan
      vb_origin = .TRUE.
    elseif ((target_origin==5).and.is_vb(5)) then 
      inten_ghost = rivb_b(i_ray2,m_x,m_z)  ! Bottom plan
      vb_origin = .TRUE.
    elseif ((target_origin==6).and.is_vb(6)) then 
      inten_ghost = rivb_t(i_ray2,m_x,m_z)  ! Top plan
      vb_origin = .TRUE.
    else
      vb_origin = .FALSE.
    endif

!----------------------------------------------------------------------------------------
! Begin integration - Compute boundary condition in Intensity 

    if (vb_origin) then  ! Virtual Boundary origin 
      inten0 = inten_ghost   
    else                 ! Open Boundary origin
! get value of j_point from the knowledge of m_x, m_y, m_z
      j_point_wall = Get_Points_Index(m_x,m_y,m_z,target_origin)

      if (j_point_wall > n_points) write(*,*) edge

      inten0       = emission(j_point_wall)  !   black body intensity      
    endif

! FOR DEBUGGING ONLY. REMOVE AFTER
    Point_origin%Temp_origin = inten0
! END

!----------------------------------------------------------------------------------------
! Loop over all the cell crossed by the ray. Starts from point at Origin
    PATHS: do i_path = n_path_ray, 1, -1 

      ! Save the indexes of the cell crossed
      Cell_path => ray(i_ray,j_point)%Passed_CV(i_path)%cub

      m_x = Cell_path%i_x
      m_y = Cell_path%i_y
      m_z = Cell_path%i_z

      surf_area = Cell_path%Surface

      inten_b = Cell_path%steradiancy       ! black body intensity in W/m^2/str

      tau_dir = exp(-Cell_path%C_abs*ray(i_ray,j_point)%length(i_path)) ! Compute transmissivity of the 
                                                                        ! path crossing the cell

      inten1  = inten0*tau_dir + inten_b*(1.0-tau_dir)
 
      ! Radiative source in the CV (net deposit of energy: minus radiative flux divergence)
      Cell_path%Source_cv = Cell_path%Source_cv+    &
                            (inten0-inten1)*cos_int(i_theta)%azimuth(i_phi)*surf_area
        
      ! Save the intensity leaving the present CV as the incident intensity 
      ! entering the next CV

      inten0 = inten1

    enddo PATHS ! end of the integration along a single ray
      
!----------------------------------------------------------------------------------------
! Save the incident irradiation along i_rayth ray

    if     ((edge==1).or.(edge==2)) then
      Array_incident(i_ray,m_x,m_y) = inten1
    elseif ((edge==3).or.(edge==4)) then
      Array_incident(i_ray,m_y,m_z) = inten1
    elseif ((edge==5).or.(edge==6)) then
      Array_incident(i_ray,m_x,m_z) = inten1
    endif

!----------------------------------------------------------------------------------------
! Compute the incident irradiation at boundary point if non virtual boundary
! If virtual boundary, set incident flux as zero

    if (.not.(vb_flag)) then 
      irrad(j_point) = irrad(j_point)+inten1*cos_int(i_theta)%azimuth(i_phi)
    else
      irrad(j_point) = 0.0
    endif
! Deallocate pointer Point_origin
    nullify(Point_origin)

  enddo Rays
        
!----------------------------------------------------------------------------------------
  return
  end subroutine int_RTE

!========================================================================================
  subroutine get_error_1D(old_irrad,new_irrad,max_variation,io)
!========================================================================================
! subroutine get_error_1D returns the maximum variation of irradiation on a real boundary
! This subroutine is used to compute the error associated with each iteration
! Update vector old_irrad
! 
! Assumed shape array. This subroutine must be set in a module.
!
! Variable passed in:
! old_irrad : (REAL) dimension(:)
! new_irrad : (REAL) dimension(:)
! io        : integer, optional. ID of the io unit
!
! Variable passed out:
! max_variation : (REAL), scalar, max_variation = max(abs(old_irrad-new_irrad)/new_irrad)
! old_irrad     : (REAL) dimension(:)
!----------------------------------------------------------------------------------------

  implicit none

! Variable passed in
  real, dimension(:), intent(in) :: new_irrad

  integer, optional              :: io

! Variable passed in and out
  real, dimension(:), intent(inout) :: old_irrad

! Variables passed out
  real, intent(out) :: max_variation
  
! Local variables
  real, allocatable, dimension(:) :: variations  

  integer, dimension(1) :: n_bound

! Get numbers of elements of irrad  

  n_bound = shape(old_irrad)

! Ensure new_irrad and old_irrad have same number of elements.
! If not call Error_MSG and shut down

  if (all(n_bound /= shape(new_irrad))) then
    if (present(io))        call Error_MSG(3,io)
    if (.not.(present(io))) call Error_MSG(3)
  endif

  allocate(variations(n_bound(1)))
  variations = 0.0
     
  where (new_irrad /= 0.0)
    variations = abs(old_irrad-new_irrad)/new_irrad
  elsewhere
    variations = abs(old_irrad-new_irrad)
  end where

  old_irrad = new_irrad

  max_variation = maxval(variations)

  if (allocated(variations)) deallocate(variations)

!----------------------------------------------------------------------------------------
  return
  end subroutine get_error_1D

!========================================================================================
  subroutine get_error_3D(old_irrad,new_irrad,max_variation,io)
!========================================================================================
! subroutine get_error_3D returns the maximum variation of irradiation on a real boundary
! This subroutine is used to compute the error associated with each iteration
! Update vector old_irrad
! 
! Assumed shape array. This subroutine must be set in a module.
!
! Variable passed in:
! old_irrad : (REAL) dimension(:,:,:)
! new_irrad : (REAL) dimension(:,:,:)
! io        : integer, optional. ID of the io unit
!
! Variable passed out:
! max_variation : (REAL), scalar, max_variation = max(abs(old_irrad-new_irrad)/new_irrad)
! old_irrad     : (REAL) dimension(:)
!----------------------------------------------------------------------------------------

  implicit none

! Variable passed in
  real, dimension(:,:,:), intent(in) :: new_irrad

  integer, optional              :: io

! Variable passed in and out
  real, dimension(:,:,:), intent(inout) :: old_irrad

! Variables passed out
  real, intent(out) :: max_variation
  
! Local variables
  real, allocatable, dimension(:,:,:) :: variations  

  integer, dimension(3) :: n_bound

! Get numbers of elements of irrad  

  n_bound = shape(old_irrad)

! Ensure new_irrad and old_irrad have same number of elements.
! If not call Error_MSG and shut down

  if (all(n_bound /= shape(new_irrad))) then
    if (present(io))        call Error_MSG(4,io)
    if (.not.(present(io))) call Error_MSG(4)
  endif

  allocate(variations(n_bound(1),n_bound(2),n_bound(3)))
  variations = 0.0
     
  where (new_irrad /= 0.0)
    variations = abs(old_irrad-new_irrad)/new_irrad
  elsewhere
    variations = abs(old_irrad-new_irrad)
  end where

  old_irrad = new_irrad

  max_variation = maxval(variations)

  if (allocated(variations)) deallocate(variations)

!----------------------------------------------------------------------------------------
  return
  end subroutine get_error_3D


!========================================================================================
  subroutine Error_MSG(arg,io)
!========================================================================================
! Thiis program is private. Prints message error depending on the value of 
! arg.
! 
! Variables passed in:
!---------------------
! arg : integer. Indicate nature of error
! io  : optional, integer. Indicates the unit for write operation. If NOT(present(io)) 
!                          then writes on the screen
!
!----------------------------------------------------------------------------------------

  use radiation_m, only : allocate_radiation_arrays

  implicit none

! Variables passed in

  integer :: arg
  integer, optional :: io

! Local variables

  integer :: DISP  

  character(100) :: MSG

  if (arg==1) MSG = 'ERROR!! Shapes of 1st and 2nd variables do not conform.'
  if (arg==2) MSG = 'ERROR!! Shapes of 1st and 3rd variables do not conform.'
  if (arg==3) MSG = 'ERROR in get_error_1D !! Shapes of 1st and 2nd variables do not conform.'
  if (arg==4) MSG = 'ERROR in get_error_3D !! Shapes of 1st and 2nd variables do not conform.'
  if (present(io)) then 
    DISP = io
  else
    DISP = 6
  endif

  write(DISP,*) trim(MSG)  

  call allocate_radiation_arrays(-1)

  stop 
!----------------------------------------------------------------------------------------
  end subroutine Error_MSG

  end module dtm_m
